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Stellingen 


behorende bij het proefschrift 


Computer Algebra Applications 


in. Econometrics 


van Arjen Merckens 





1. 


Slechts enige basiskennis ‘is nodig om in principe in staat te zijn 


computeralgebra te bedrijven in een hogere programmeertaal. 


De wijze van evaluatie van de rangconditie voor simultane vergelijkingen 
met meetfouten zoals gegeven door Hsiao in het Handbook of Econometrics is 
onjuist en kan tot verkeerde conclusies leiden met betrekking tot de 
identificatie-situatie van de parameters van een model. ({Merckens, A. en 
P.A. Bekker, 1988; tevens Hoofstuk 3 van dit proefschrift] 


. Definieer een parameter—matrix/vector als een matrix/vector waarvan elk 


element een afzonderlijke parameter representeert. Als op de elementen van 
een parameter-matrix F en parameter-vector f slechts exclusie—restricties 
rusten, en F een matrix van volledige kolommenrang is, dan is f lineair 
afhankelijk van de kolommen van F d.e.s.d.a. er n>0 te vinden is, zodat er 
een submatrix F,, van n verschillende kolommen van F bestaat waarvoor geldt 
dat. exact n rijen van (F,, f) niet gelijk aan een nul-rij zijn. [Zie 
Theorem 5.3.1 van dit proefschrift] 


. De "repeated conditioning” methode kan gebruikt worden om de verwachting 


van polynomiale functies van normaal verdeelde matrices analytisch te 
berekenen. [Zie Hoofdstuk 7 van dit proefschrift] 


. Twee modellen zijn equivalent als zij dezelfde set van kansverdelingen 


voor de waarneembare vector impliceren. In dat geval kunnen de data nooit 
voldoende informatie verschaffen om onderscheid tussen de twee modellen te 


maken. Beschouw nu een systeem van simultane vergelijkingen 
By = [x + €, 


met B een non-singuliere mxm—matrix, [ een mxk-matrix en de storingen € 
normaal verdeeld en onafhankelijk van x. Als op de parameter—matrices B en 
I slechts exclusie-restricties rusten, dan kan de procedure voor het 
bepalen van een minimale parameterisatie van BT tevens gebruikt worden 


om de equivalentie van dergelijke modellen te onderzoeken. 





10. 


i. 


. De grote bekendheid van het verschijnsel computervirus heeft als gunstig 


neven-effect dat men zich de kwetsbaarheid van computers beter realiseert 
en men meer tijd besteedt aan het maken van backups, zodat eventuele 
"hardware failures” minder verantwoordelijk zijn voor het verloren gaan 


van grote hoeveelheden werk. 


. Een eenvoudige en effectieve maatregel om te hard rijdende automobilisten 


. aan te pakken is het direct in beslag nemen van het rijbewijs voor het 


aantal dagen dat er kilometers per uur te hard gereden is. Het 
afschrik—-effect voor het kwijtraken van hun mobiliteit is dan zo groot dat 


de pakkans even laag kan blijven als nu het geval is. 


. Het is financieel gezien attractief om een "Ph.D. student” in Nederland te 


zijn. 


Het getuigt van immoreel gedrag om Shareware Software te zien als 


Freeware. 


Thuiswerken moet waar mogelijk gestimuleerd worden: het is een misvatting 
dat thuiswerkers over het algemeen niet produktief zijn, zodat er, bijv. 
door de vermindering van het aantal files en de vermindering van de 


reistijd, grote economische winsten zijn te boeken. 


”Wanneer der wetenschap weet wat der toekomst brengen gaat, dan kan der 
wetenschap niet verder gaan. Het is niet wetenschappelijk om vooruit te 
kijken” ~ Prof. Prlwytzkovski in De Niks uit de bundel Prawl!! Der 


Hemeldonderweder van Marten Toonder. 





CHAPTER 1 


INTRODUCTION 


In this study we will introduce computer algebra in econometrics and 
statistics. In short, computer algebra programs do symbolic mathematics: 
computer algebra programs are not only able to calculate numerical expressions 
without rounding error, but also able to manipulate formulas without numerical 
values; any algebraic object can be represented exactly in the memory of the 
computer. , 

Computer algebra is not yet commonly used in solving econometric problems, 
though it could be a useful tool for applied as well as for theoretical 
econometricians. An applied econometrician could, for example, let a computer 
algebra system — an automatic system for symbol manipulation — compute the 
derivative of some object function and even let the resulting formula 
automatically be programmed - error free. A theoretical econometrician could 
use a computer algebra system as a tool for performing complicated or tiresome 
manipulations on some analytical expressions. 

The aim of this study is to investigate the integration of computer algebra 
in econometrics. In other words, we will investigate if computer algebra is a 
useful tool in tackling a diversity of econometric and statistical problems. 
In each chapter, except chapter 2, we will give the econometric or statistical 
context of a problem and give algorithms to solve the problem. under 
consideration. All these algorithms have in common that they have to be 
implemented in computer algebra programs. 

In this thesis hardly any attention is given to existing computer algebra 
systems; all problems in this thesis have been handled by special purpose 
computer algebra programs written in the higher-level programming language 
Pascal. An advantage of writing computer algebra applications. in Pascal is 
that it yields better insight in the limitations of computer algebra, since 
the algorithms used are not black boxes, but self implemented algorithms. 
Another advantage is that it not restrictive: every manipulation or every 
algorithm can be programmed, whereas in computer algebra systems an algorithm 
needed might not be available, or a manipulation might not be defined. An 
advantage of computer algebra systems over a higher-level programming language 
is, that the user does not have to know how symbols are stored in the memory 
of the computer. Also, computer algebra systems contain standard libraries of 
algorithms for manipulation of expressions. 


Chapter 2 will discuss how computer algebra programs can be programmed in a 
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higher-level programming language like Pascal. It will do so by discussing all 
the algorithms needed to build a basic, simple computer algebra system. 

Chapters 3, 4 and 6 handle the identification problem of several models. In 
all these chapters it is assumed that the observations follow the normal 
distribution and are identically and independently distributed. So the 
identification problem is equal to the question whether the parameter values 
can be uniquely determined from the first- and second—order moment equations. 
Usually rank conditions are given for the identification of the parameters of 
structural models. Often the applied econometrician is left with the question 
how to actually evaluate these rank conditions. In chapter 3 a computer 
program (ERA) will be presented that computes the rank - and null—space — in 
these conditions for a variety of models. 

Hsiao (1983) and Geraci (1976) presented rank conditions to characterize 
the identification problem of simultaneous equation models with measurement. 
error. It will be shown in chapter 3 that these rank conditions are not 
suitable for evaluation: they might not yield the correct result. An 
alternative rank condition is derived and a computér program (ERASIM) is 
presented that uses the parameter restrictions as input in order to 
characterize the identification situation of the individual parameters in the 
output. 

Chapter 4 investigates the identification of restricted factor analysis 
models. In this chapter also a general rank condition is derived and a 
computer algebra program is presented (ERAFAC) that can automatically 
construct and evaluate the Jacobian matrix used to characterize the 
identification situation of the individual parameters. 

Bekker and Dijkstra (1990) discuss the precise number of restrictions on 
the reduced form of a simultaneous equation system as implied by exclusion 
restrictions in the structural form. In chapter 5 an algorithm, to a large 
extent due to Bekker, and a computer program (MINPARAM) are presented to 
characterize these restrictions with a minimal parameterization of the reduced 
form, ie. a parameterization taking into account ‘the restrictions on the 
reduced form and using as few as possible parameters. 

A minimal parameterization of a matrix like B"'I is also used in chapter 6. 
This chapter handles the identification problem in general structural linear 
models, more commonly known as LISREL models (Jéreskog, 1984). A number of 
rank conditions are derived and also a number of computer algebra programs are 
presented able to construct and evaluate these rank conditions. The first part 


of this chapter deals with the identification problem of specific submodels of 
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the general LISREL model, whereas the last part of this chapter handles the 
identification problem of LISREL models with general linear restrictions on 
the parameter matrices. 

In chapter 7 the theoretical problem of the expectation of higher-degree 
functions of normally distributed matrices is investigated, an area where 
there is a certain amount of recent literature, see for example Neudecker and 
Wansbeek (1987). It is shown that the expectation of these functions can be 
evaluated algebraically by a method that is called repeated conditioning. A 
computer program is presented that can handle a large diversity of this type 
of problems. 

Chapter 8 introduces some sort of computerized scratchpad, capable of 
performing matrix algebra in the sense that symbols denote matrices and the 
matrices cannot be filled with specific elements. Matrices can have certain 
characteristics, like idempotency, or a matrix can be substituted for a matrix 
expression. This scratchpad ~ MATRAL - can be used to compute a wide variety 
of symbolic matrix expressions. Here it is applied to compute the variance of 
quadratic unbiased estimators of the error—-component model with incomplete 
panels, as discussed by Wansbeek and Kapteyn (1989). 

The seven main chapters that this thesis embodies will show that computer 
algebra is a useful instrument in solving certain theoretical and practical 
problems in econometrics. Of course, the application of computer algebra is 
not limited to the scope of problems handled by this thesis. 








CHAPTER 2 


COMPUTER ALGEBRA 


2.1 Introduction 


Computer algebra can be seen as that part of computer science which designs, 
arialyses, implements and applies algebraic algorithms. It does not only deal 
with non-numerical computing: any algebraic object can be represented exactly 
in the memory of a computer, so that algebraic computations can be performed 
without loss of precision and significance. 

A long history of formula manipulation by computer exists in several 
fields, like physics. In the early sixties the first special purpose computer 
algebra programs were constructed. Complete systems to perform computer 
algebra like Reduce (1967) and Macsyma (1973) arrived. These systems were 
capable to perform high school mathematics: symbolic differentiation, Taylor 
series, and even integration became possible at the push of a button. In 1979 
Stoutemeyer states ‘Symbolic computation comes of age’. However, in fields 
like econometrics computer algebra was hardly heard of, and even at this time 
it is barely exploited. However, with the arrival of more powerful computers 
and because more people are getting aware of the existence of computer algebra 
systems, this is changing rapidly. 

Computer algebra systems are automatic systems for symbol manipulation. 
These systems enable a user not only to solve standard problems, but also to 
implement algebraic objects in a computer program without much effort. A brief 
discussion with some examples will be given in section 2.2. 

Sometimes these computer algebra systems are not at hand or are not capable 
to tackle the problem in mind..The aim of this chapter is to show how computer 
algebra problems can be tackled by using a higher-level programming language 
like Pascal. In order to be capable of manipulating a formula in Pascal, some 
basic knowledge about pointers, trees and linked lists is necessary. Also some 
basic algorithms are needed. Section 2.3 to 2.10 will (briefly) discuss these 


topics. In section 2.11 the information in the previous sections will be 


combined, which yields as a result a Very Simple Computer Algebra System 
(VSCAS). Section 2.12 concludes. 
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2.2 Computer Algebra Systems 


Using a computer algebra system (CAS) will suffice for solving many types of 
problems. Analytically solving systems of linear equations or computing the 
derivative of a function is often straightforward using a CAS; these are built 
as an integral part of a CAS. Most systems also include a programming language 
which can be used to write computer programs that are able to store and 
manipulate algebraic expressions. Examples of some well-known systems are 
Reduce, Macsyma, Maple and - rather new - Mathematica. The last computer 
algebra system is even discussed in an economic journal by Belsley (1989). 

The output a CAS gives may sometimes surprise a novice user: when a user 
types (a+5)*(b+2) he may get the answer in expanded form: axb + bxx2+2xa+2*b, 
in factored form: (a+b)*(b+2) or in recursive form: a%(b+2) + biek2 + 2xb. However, 
often he can specify which form he prefers. 

In this thesis it is not intended to give a complete overview of the 
capabilities of some (or all) computer algebra systems. By giving a few 
self—-explaining examples some notion is given of the easiness to use and the 
power of CAS’s. The CAS used for the examples is Reduce. For a thorough 
description on most CAS’s, confer van Hulzen and Calmet (1982), Char et al. 
(1983, 1985) and Wolfram (1988). 


Example 2.2.1 Solving equations. 


INPUT: 
equation] := 3axx+4y-2z2 =k 
equation2 := 7x-4axy+4z = 1 
equation3 := 3x+6xaxy+5z = m 


SOLVE({equation1, equation2, equation3},{x,y,z}) 


OUTPUT: 
{ xa_llxAxkK + 3xAxL + 2%AxM + 5eL -— 4%M 
Fa ft RAL = lp ea la ; 
334A? + 27KA + 23 


ye _ _L5#AwL - 12%AWM - 234K + GL - 144M 


Zan-280_Ko_t_2EO AD = ARSE <2 ht _a ake } 
2x(33KA7 + 27%A + 23) 
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Example 2.2.2 Integrating a function. 


INPUT: 
int(4*xxlog(x)+xxcos(xX)—x*x«siN(X),X) ; 


OUTPUT: 
— (SIN(X)*X — COS(X)*X” + COS(X) — 2xLOG(X)X? + X?) 


Example 2.2.3 Computing the Taylor series. 


INPUT: 
(Note: SUB means substitute, DF means derivative) 
Algebraic procedure Taylor(functionx, x, point, degreen) ; 
SUB(x=point, functionx) + 
FOR k := l:degreen 
SUM( SUB(x=point, DF(functionx, x, k)) 
*(x—point)xeKk / 
(FOR j := 1:k product j) 
i 


tl := Taylor(Ex«x,x,0,4) ; 
ON DIV,RAT 
t2 := Taylor(In(1+y),y,0,6) 


OUTPUT: 


te Xt 44x? + 124X? + 244K + 24 


t2=2-(—3-#¥°— FWY" a oe ais —5-#¥*-Y) 
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2.3 Computer Algebra in Pascal . 


In order to be able to understand how computer algebra programs can be 
constructed in a conventional programming language like Pascal, basic 
knowledge about pointers and data structures, such as linked lists and trees, 
is necessary. The next three sections will be dedicated to this. In the 
section on trees, a data structure for manipulating algebraic expressions will 
be given that will be used in the rest of this chapter. Section 2.7 will give 
algorithms for calculating without loss of precision. Section 2.8 will give 
algorithms for addition and multiplication of algebraic expressions. In order 
to represent expressions in a unique way, a sorting method has to be used. A 
suitable sorting algorithm will be described in section 2.9. Before 
expressions can be manipulated, these expressions have to be read and stored 
in memory first. However, this is more complex than reading and storing an 
integer. Section 2.10 will show a parser that can handle this. 

The aim of this chapter is the construction of a very simple CAS. At the 
end of this section an outline of a very simple CAS programmed in Pascal will 
be presented, using data structures and algorithms as given in this chapter. 
Computer algebra will lose much of its magic — as some people see it — when it 
is clear that it is a mere combination of data structures, a parser, a sorting 
method and a few algorithms. For extensive descriptions on some of these 
topics see for example Aho, Hopcroft and Ullman (1983), Knuth (1981) or 
Sedgewick (1988). 


2.4 Data structures 


In order to avoid any possible misunderstanding about the concept of data 
type, abstract data type and data structure, the definitions of these terms 
will be given. 


Definition 2.4.1 
The data type (or just type) of a variable is the set of values a variable 
may assume (the domain). 


Definition 2.4.2 
An abstract data type (ADT) is a mathematical model with various operations 
defined on the model. 
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Definition 2.4.3 
A data structure is a collection of variables, not necessarily of the same 
type, connected in various ways. Data structures are used to represent a 


mathematical model of an ADT. 
Example 2.4.1 


- The set of values (TRUE, FALSE) is the data type of a boolean variable. 

- A set of vectors together with the operations addition, subtraction, 
transposition and multiplication is an example of an ADT. 

- An array is the simplest form of a data structure; a one-dimensional array 


can be seen as a representation of a vector. 


A cell is the basic building block of a data structure. A cell can be thought 
of as a box, capable of holding values drawn from some data type. Another 
common mechanism for grouping cells in a programming language is what is 
called in Pascal the record structure. A record is a cell that is made up of a 
collection of cells, called fields of possibly dissimilar types. 


Example 2.4.2 A record structure. 


type 
recordtype =record 
number: integer; 
character : char; 


end; 


Most higher-level programming languages are also capable of representing 
relationships between cells using pointers. A pointer is a cell whose value 
indicates another cell. In pictures of data structures an arrow going from 
cell A to cell B will indicate that A is a pointer and points to B. In Pascal 
a pointer—-variable to a cell of type recordtype can be created by the 


declaration: 


var 


ptr: “recordtype; 


| 
I 
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The expression ptr* denotes the value (of type recordtype) in the cell pointed 
to by ptr. Diagram 2.4.1 shows the situation after the statements 


new(ptr); ptr’.character := A; ptr°.number := 3; 


Diagram 2.4.1 
SHEE 


The statement new(ptr) allocates an amount of memory capable of holding the 
information of type recordtype; a statement dispose(ptr) will free the 
allocated memory for other use. This type of memory usage is called dynamic 
memory allocation, in contrast to static memory allocation, which is memory 


reserved for variables at compilation time. 


2.5 Linked lists 


Mathematically, a list is a sequence of zero or more elements of a given type. 
Often such a list is represented by a sequence of elements, separated by 


comma’s: 
Gy, Ga, .. 1 An n > 0, each a; of the same cell type. 


By defining a set of operations like insert, delete and locate on objects of 
type list, an abstract data type will be formed. It is easy to see that the 
data structure array can be used to represent lists. 

Another data structure able to represent a list is the so called linked 
list. A linked list is the most basic form of a dynamic data structure. It is 
a sequence of nodes in which each node is connected (linked) to the node 
preceding it. Each node contains data and a pointer to the next element. See 


diagram 2.5.1 for an example. 


Diagram 2.5.1 


Tg AE i es Ce 


head 
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A pointer value of [J denotes the end of the list, or, in Pascal nil (points to 
nothing); head denotes the start of the list. 


The basic cell of a linked list can be defined by the definition 


type 
nodeptr = “node; 
node = record 
number: integer; 
character : char; 
next : nodeptr; 
end; 


Assuming that a list should be kept sorted, one of the advantages of linked 
lists over the array implementation of list is that insertion and deletion of 
a cell is fast, as will be shown by example 2.5.1. A disadvantage is that 
direct access to a cell is not possible: to find the value of cell k in the 
list, K-1 pointers have to be followed. 


Example 2.5.1 Array insertion versus linked list insertion. 

Diagram 2.5.2a shows the array insertion of a cell, whereas diagram 2.5.2b 
shows the insertion of a cell in a linked list. 

The array implementation gives direct access to cells, so a method like 


bisection can be used in order to find the location where the cell |’7” 


should be inserted. In the linked list implementation no direct access to 
cells exists, so the correct location to insert a cell can only be found by 
sequential or linear search. So search cost is more expensive in the linked 
list implementation. 

In the array implementation n-k cells have to be copied, whereas in the 
linked list implementation only 2 pointers have to be (re)set. 

On the average the array implementation has *log(n) comparisons, and 1+n/2 
cell assignments. The linked list implementation has n/2 comparisons, 1 cell 
assignment and 2 pointer assignments. Considering that cell assignments are 
more expensive than (field—) comparisons, it is clear that insertion using the 


linked list implementation is less expensive than the array implementation. 


sal 
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Diagram 2.5.2a. Array insertion of cell [7s |. 


bed 


A 
= 


A 


a 
ms 
Le | 
=" 
> 
~ 
| 
~ 


rod 


to 
to 


Ww k+1)’T’| 3 


tow 

re 

_ 

Ee 

tow 

+ Co 

a 
> 


k+2|U? k42|U" 


\ 
El 
ie 
el 
i 


n+1]|’?Z’ n+1|'?2’ 


I. The original situation: the cell must be inserted between cell 


k and k+1 (after ’R’, before ’U’). 
II. Copy cell n to cell n+1, n-1 to n, ... K+1 to k+2. 


Ill. Store value |’7” in cell k+1. This completes the insertion. 
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Diagram 2.5.2b Linked list insertion of cell [7s]. 


HEN alc gy ok ald 
weal el Gal ood OE cn 





a, pee eee ae a) Pela ela 





I. The original situation: the cell must be inserted behind cell 
[@] 9] |}. This k-th cell can be found by following k-1 pointers. 

II. Create a new cell, store the value in it and let it point to 
ceo] 2 | |. 


III. Let cell k point to the new cell. This completes the insertion. 


13 
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2.6 Trees 


A tree imposes ‘a hierarchical structure on a collection of items. Familiar 


examples of trees are genealogies and organization charts. Trees are common in 


many areas of Computer Science; for example, trees are used to organize 
information in database systems. Definitions 2.6.1 - 2.6.4 will define some 


terms commonly used when talking of trees. 


Definition 2.6.1 
A vertex is a simple object (also called node) that can have a name and can 


carry other associated information. In terms of the previous section: a cell. 


Definition 2.6.2 
An edge is a connection between two vertices. In terms of the previous 


section: a pointer. 


Definition 2.6.3 
_A path in a tree is a list of distinct vertices in which successive 
vertices are connected by edges in the tree. 


Definition 2.6.4 

A tree is a non-empty collection of vertices and edges. One vertex in the 
tree is designated as root. A tree must have exactly one path between the root 
and each of the other vertices in the tree. 


As is clear from these dbfinitions, a (linked) list is a simple example of a 
tree. 

Trees may be used for representing expressions in memory. For example, the 
arithmetic expression (a+b)*x(a—b) could be represented in an expression tree 
as drawn in diagram 2.6.1. : 

When manipulating formulas by computer, it is mandatory that equivalent 
expressions are recognized as such. Unfortunately, equivalent expressions may 
have different expression trees which complicates the recognition of equal 
expressions. Consider for example diagrams 2.6.1 and 2.6.2. 
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Diagram 2.6.1 Expression tree for (a+b)x(a—b). 


root 


interior 
nodes 





Diagram 2.6.2 Expression tree for (axa) —(bxb). 


Ugo ae oe OG 


We all know that the expressions (a+b)x(a—b) and (axa)-(bxb) are equivalent; a 
computer program must first normalize these expressions before the equality of 
these expressions can be recognized. This means that all equivalent 
expressions should be represented by the same expression tree. An obvious 
approach to normalize expressions is to remove the influence of parentheses, 
sort factors, sort terms, add equal terms and eliminate zero terms. A suitable 
algorithm can be found in de Jong (1989, p.69-74). 

Expression trees are just one way of storing algebraic expressions: many 
other data structures can be defined. In diagram 2.6.3 an example is given of 
a data structure than can only store expressions in expanded form. The 
operations * and + are implicitly defined within the structure: an arrow 
(pointer) in a horizontal direction denotes an addition and in vertical 


direction a multiplication. 
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Diagram 2.6.3 Representation of a?—b?. 


Beclis 
a 


In Pascal records and multiple pointers are used to define the data structure 
as indicated in diagram 2.6.1. See example 2.6.1 for a formal declaration of 





this data structure. 


Example 2.6.1 A data structure definition in Pascal 


type 
factorptr = “factortype; 
factortype = record 
character: char; 
next: factorptr; 
end; 
expressionptr = “expressiontype; 
expressiontype = record ; 
scalar : integer; 
factor : factorptr; 
next : expressionptr; 


end; 


When this data structure is used and input should be allowed to be of a more 
general form than the expanded form (e.g. input expressions may contain 
parentheses), then at parse time — i.e. when reading an expression from input 
- the problem of storing expressions that contain parentheses has to be 


solved. This will be discussed .in the section on parsing. 
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2.7 Arithmetic 


When using computer algebra systems, the multiplication of two numbers will 
(and should) have an exact outcome. However, in conventional programming 
languages, like Pascal, precision is limited. So, when computer algebra has to 
be performed using a conventional language, some algorithms have to be 
implemented that can add, subtract, multiply and divide numbers that have 
arbitrarily high precision. For simplicity, assume we are working with 
integers instead of with numbers with an imbedded radix point. Because of the 
importance of these algorithms in computer algebra a rather lengthy 


description will be given below. The algorithms we will discuss are 


1. addition or subtraction of n—place integers giving an n-place number and 
a carry. 

2. multiplication of an n-place integer with an m-place integer giving an 
(m+n)—place answer. 

3. division of an (m+n)-place integer by an n-place integer giving an 
(m+1)-place quotient and an n-place remainder. 


The algorithms will do operations 1, 2 and 3 in radix—b notation, where 5b is 


any given integer > 2; so, 


(0.0@p-34n-2%n14n)p = ---+Op_gb” +a,-2b7+ an.b'+ a,b°, and 
(520), = 5.6°+ 2.6'+ 0.6° = 192. 


In all these algorithms we shall assume that we deal with non-negative 


numbers. 
Algorithm 2.7.1 Addition of n—-place integers [from Knuth (1981, p.251)]. 


Given non-negative n-—place integers (t,U2..Un)y and (0,V2..U,),, this 
algorithm forms their radix—b sum, (wpw,W2..W,)p, With wy the carry which will 
always be equal to 0 or 1. 


1. [Initialize] Set j <n, k <0. 
[Add digits] Set w; <— (u;+v;+k) mod b, and k <- (u;+v,;+k) div b. 
{Loop on j ] Decrease 7 by one. If j > 0, go back to step 2; otherwise set 
Wg <— k and terminate the algorithm. 


17 


Computer Algebra 


An algorithm for subtracting an n-place integer from another is similar to 
addition and the corresponding algorithm can be found in Knuth (1981). These 
two algorithms are a representation of conventional arithmetic: in principle, 
when using paper and pencil, the same manner of calculating is used. An 
example below shows for the addition of two numbers using algorithm 2.7.1, 


with the values of j, k, u,;, vj and w, given in each step. 


Example 2.7.1 Addition of (134), with (424)s. 


Step j Us; v; k Wyo Wy Wz W3 
> & d@ dea 3 3 
2 2 3 2 421 1 3 

2 1 1 4 1 1 1 3 
3 0 1 1 1 1 3 


So, the summation of these two numbers is w = (1113). 


Algorithm 2.7.2 Multiplication of an n—place integer with an m-—place integer 
[from Knuth (1981, p.253)]. 


Given non-negative integers (tyta..ty), and (0,V2..Um),, this algorithm 
forms their radix—b product, (wW,W2.-Wmin)p Addition of the partial products 


will be.done concurrently with the multiplication. 


[Initialize] Set Wmsi)Wms2)-->Wmen all to zero. Set j <— m. 

[Zero multiplier?] If v; = 0, set w; <- 0 and go to step 6. 

[Initialize i] Set i <n, k < 0. 

[Multiply and add] Set ¢ <— u;xv; + w,,; + k; then set w,,; <— ¢ mod b and 
k <— |t/b|. 

5. [Loop on #] Decrease i by one. Now if i > 0, go back to step. 4; otherwise 


a YP 


set w; <— k, 
6. [Loop on j] Decrease j by one. Now if j > 0, go back to step 2; otherwise 


the algorithm terminates. 
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Example 2.7.2 Multiplication of 914 by 84, assuming that radix—b = 10. 


Step ij uv; t k WW, Wye Wy Wy Ws 
5 3. 2 4 4 161 0 O 6 
5 2 2 1 4 05 0 0 5 6 
5 cd -2 9 4 36 3 6 5 6 
6 0 2 4 36 3 3 6 5 6 
5 3 41 4 8 37 3 3 6 7 6 
5 2 1 1 8 i7 1 3.7 #7 6 
5 1 #1 9 8 6 7 - 6 7 7 6 
6 0 1 8 76 7 7 6 7 7 6 


There exist faster methods for the multiplication of two large integers. 
However, this method has the advantage of simplicity. 


The third and most difficult algorithm is the long division. When dividing 
two integers using paper and pencil, a certain amount of guesswork is 
involved. Before getting a functional algorithm this guesswork has to be 
explained. Long division breaks down in simpler steps each of which is the 
division of an (n+1)-place number u by the n-place divisor v. The remainder r 
after each step is less than v, so we may use rb+(next place of dividend) as 
the new u in the succeeding step. For example, if we have to divide 4183 by 
62, then we first divide 418 by 62 getting 6 and an remainder of 46. Then we 
divide 463 by 62 getting 7 and a remainder of 19. So we have a quotient of 67 
and a remainder of 19. As a consequence an algorithm for long division can be 


constructed if we can construct an algorithm to solve this problem: 


Let u = (UpUjUg..U,), and v = (v,05..U,), be non-negative integers such 
that u/v < b. Determine q = |u/v]. 


An obvious approach to solve this problem is to make a guess about qg based on 
the most significant digits of n and v. 


Upb+U, 
¢@ = min , b-1]. 
vy 


It can be proved that this g is never too small and in particular: 
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Theorem 2.7.1 
If v, > [b/2], then 9-2 < q < @. (Proof: see Knuth (1981), p.257) 


So @ is a good estimate if v, > |b/2]. This condition can ‘be seen as a 
normalization requirement. This normalization can be done by multiplying both 
u and v by |b/(v,+1)|; this does not change the value |u/v|. Now we are ready 
to present the algorithm for long division. 


Algorithm 2.7.3 Division of am m+n- place integer by an n-place integer [from 
Knuth (1981, p. 257,258)}. 


Given non-negative integers (u,U.-Umin)p aNd (v,V2..0,), Where v, # 0 and 
n > 1, this algorithm forms |u/v| = (¢o91--¢m)y and the remainder u mod v = 


(Tyo. -Tr)p 


1, [Normalize] Set d <— [b/(v,+1)|. Then set (uto..Umin)y equal to 
(Ujtg-.Umin)p times d, and set (v,v2..¥,), equal to (v,V2..U,), times d. 
[Initialize j] Set 7 < 0. 

[Calculate @] If uj;=v,, set @ <— 5-1; otherwise set 9 <— [(ujb+t6541)/v4] - 
Now test if upg > (ujb+u;41—9v,)b+u;42; if so, decrease 9 by 1 and repeat 
this test. ; 

4, [Multiply and subtract] Replace (UjUj41--Uyan)n bY (UjUj41--Ujan), minus G 
times (v,02..U;)p. 

5. [Test remainder] Set 9; <— q. If the result of step 4 was negative, go to 
step 6; otherwise go to step 7. 

6. [Add back] Decrease q; by 1, and add (0v,v9..v,), to (Ujtej41-- Ujan)D- 

[Loop on j] Increase j by one. Now if j < m, go back to 3. 
[Unnormalize] Now (9091--%m)» is the desired quotient, and the desired 


remainder may be obtained by dividing (tmtm41--Umin)p by d. 


Since this algorithm is rather complicated, a thorough look at example 
2.7.3 might help to understand the algorithm. \ 
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Example 2.7.3 Long division of (4234), by (134)s. 


Step 1,2. d = [5/2] = 2, j=0, u = (upuyugugts)s = 2%(04234), = (14023),, 
DV = (V,U203)5 = 2%(134)5; = (323)5. 
Step 3. @ = [9/3] = 3. 


Veg = 6 > (uUgb+u,—Guv,)b+u, = 0, so J = 3-1 = 2. 
Step 4. (tpt ttgttg); = (1402),—2%(323), = (201)5. 
Step 5,7. gg = 0, j = 1, and u = (uyugugt,), = (2013)s. 
Step 3. G = [10/3] =-3. 
UG = 6 < (uyb+U2—Gr,)b+us = 6. 
Step 4. (tu Uttgttg)s = (2013),—3%(323), = (2013);-(2024), = —(0011)s. 
Step 5. q, = 3. 
Step 6. gq, = 3-1 = 2. (uytgtguy); = —(0011), + (0323), = (0312),. 
Step7. j=2>me=1. 
Step 8.  (22)s is the desired quotient. The remainder is obtained by dividing 
(0312); by 2 and equals (131)s. 


2.8 Manipulation of expressions 


In the previous section we discussed algorithms for arithmetic on numbers, 
in this section we will discuss algorithms for addition, subtraction and 
multiplication on expressions. We will use the type definition of 


expressiontype as in example 2.6.1: 


expressiontype = record 
scalar =: integer; 
factor : factorptr; 
next : expressionptr 


end; 


Scalar should have a data type that can be manipulated with unlimited 
precision. This can be done by defining a data type to which the algorithms in 
the previous section can be applied; to simplify somewhat the algorithms to 
come, we assume that arithmetic on integers can be done with unlimited 


precision. 
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Definition 2.8.1 
xq is the i-th element in the list after element x. If the i-th element 
does not exist, then x(j) = nil. xq) = x. 


~ 


Algorithm 2.8.1 Addition, of two expressions. 


Given, two expressions x and y of data type expressionptr, this algorithm 


forms their unsorted sum z. 


1. [Initialize] Allocate memory for h; set z <— h. h is a variable which 
contains the head (start) of the list. 

2. [End loop] If. x = nil then go to step 5. 
[Add x to z] Allocate memory for 2) . Set z <— Zq)s 
2’. factor <— x*. factor , z*.scalar < x*. scalar. 

4. [Loop on x] Set x < X(1) and go to back to step 2. 
[End loop y] If y = nil then go to step 8. 
{Add y to z] Allocate memory for Za) - Set z <— zy) , 
2. factor <— y*. factor , z*.scalar < y*.scalar. 
[Loop on y] Set y <— yj) and go to back to step 5. 
[Finish] Set z < hy) and free memory taken by h. 


Example 2.8.1 will show how two expressions are added using this algorithm. 


Example 2.8.1 Addition of two expressions x and y, with 


«= ES ST ony = TS. 
Cm tr CH ODM 


3: hag = [8 7 Ti 
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3: hay = 
fm oh 


_¢ iil 
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6: hay = Bl 4-GT 4-4 TT | 


Cm Gh CHE 
Ce Tl 


6 hg) = [BT 1-41 4-1 


mm | s a | 
! 


| fal 


8 z= hay 


Subtracting expression y from expression x is almost equal to algorithm 
2.8.1, except that in step 6 the assignment z*.scalar < y*.scalar has to be 
changed in z*. scalar <— — y*. scalar. 


Algorithm 2.8.2 Multiplication of two expressions. 


Given two expressions x and y of data type expressionptr, this algorithm 
forms their unsorted product z. 


[Initialize] Allocate memory for h, set z < h and s < y. 
[End loop x] If x = nil then go to step 10. 

[Reset y] Set y < s. 

{End loop y] If y = nil then go to step 9. 

[Reserve memory] Allocate memory for Za) and set z < 2). 
[Multiply scalars] 2°. scalar <— x*.scalar * y*. scalar. 


SS Oe iS iS me 


(Multiply factors] z*. factor <— x*.factor * y*.factor; of course a 
procedure has to be written that establishes this factor multiplication. 


% 


[Loop on y] Set y <— yy) and go back to step 4. 
9. [Loop on x] Set x <— xi) and go back to step 2. 
10. [Finish] z <— hi), and free memory taken by h. 


The example below will show how two expressions are multiplied using this 
algorithm. 
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Example 2.8.2 Multiplication of two expressions x and y, with 
dea BT I~ and y = GT Sa _- 
el | (ef) Le mz | 
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All these algorithms give as result an unsorted and thus not completely 
normalized expression. In’ order to complete normalization the terms within an 
expression have to be sorted, the equal terms added and the zero terms 


deleted. £ 
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2.9 Sorting 


As was mentioned at the end of the previous section, the resulting expression 
has to be sorted in order to establish a normalized expression. There are many 
sorting methods available. One of the most efficient methods when sorting 
random elements is Quicksort. However, Quicksort can only be efficiently 
implemented if the elements to be sorted are directly accessible. Therefore, 
implementing Quicksort for lists (like our data type expressionptr) is not 
possible. Another sorting method is extremely well suited for sorting lists: 
Mergesort. Mergesort is based on the principle that two sorted lists (or: 
arrays, files) can easily be merged into one sorted list. It is also based on 
the divide-and-—conquer paradigm: divide the list in half, sort the two halves 
(recursively) and merge the two halves into one sorted list. 

Another sorting method capable of sorting lists is Bubblesort: keep passing 
through the list, exchanging adjacent elements, if necessary; when no 
exchanges are required on some pass, the list is sorted. However, Bubblesort 
is known to be one of the most inefficient sorting algorithms. There are more 
sorting methods than the three named here, but when sorting lists, Mergesort 
is known to be one of the most efficient sorting methods. 

In this thesis a non-recursive version of Mergesort will be used: scan 
through the list performing 1 by 1 merges to produce sorted sublists of size 
2, then scan through the list performing 2 by 2 merges to produce sorted 
sublists of size 4, then scan through the list performing 4 by 4 merges to 
produce sorted sublists of size 8, etc., until the list is sorted. 

In this section a non-recursive mergesort algorithm for lists will be 


given. First a merging algorithm is needed. 
Algorithm 2.9.1 Merge two lists. 


Given two ordered linked lists x and y, this algorithm merges these two 


lists into one ordered list. 


1. [Initialize] Set z <— h. It is assumed memory has been allocated to h by 
an initialization procedure. 

2. [Compare] If x*.factor < y’.factor then set 2(1)<— ‘x, z2<— x and x <— X1) 
otherwise set Z1) <- y, Z <— y and y <— a). 

3. [End of lists?] If x # nil or y # nil then go back to step 2. 
[Finish] The start of the merged list can be found in hy). 
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In step 2 the < sign indicates a procedure performing a comparison on the two 
factors named. Note that in this algorithm the original lists x and y are 
destroyed. 


This merging algorithm (2.9.1) will be used in the following "bottom-up’ 


mergesort algorithm. 
Algorithm 2.9.2 Non-recursive Mergesort. 
This algorithm sorts a randomly ordered linked list z. 


[Initialize] Allocate memory for h and set hy) <— z and n <— 1. 
{Reset r and z] Set r <— hy) and z — h. 

[Start of list a] Set t <— 7, a < t and t — ty). 

[Start of list b] Set b < ty). 

[End of list a] Set t() < nil. 

[End of list b] Set  <— b, t <— tinay, 7 <— tay, tay <— nil. 
[Merge sorted sublists a and b] 24) <— merge(a,b). 

[Skip 2n sorted elements] z <— Zan). 

(Loop] If r # nil then go back to step 3. 

[Next merging phase] n <— 2xn. 

[Loop] if a # h then go back to step 2. 


Sere A Sf wm 


RP = 
PFS 


[Finish] z <— hi) and free memory taken by h. 


Note: function merge in step 7 returns a pointer to the sorted list produced 
by merging lists a and 5 (hy) as given in algorithm 2.9.1). 

This mergesort algorithm is demonstrated in diagram 2.9.1: the first row gives 
the unsorted list. The second row shows that two sublists of length 1 (”u” and 
*n”) are merged, the third row shows that again two sublists of length 1 (”s” 
and ”o”) are merged, etc. The eight row shows that two sublists of length two 
are merged ("nu” and ”os”) etc., the eleventh row show that two sublists of 
length 4 are merged (”"nosu” and “dert”), and the last row shows that a sublist 
of length 8 ("denorstu”) and a sublist of length 4 ("ilst”) are merged, so 


that all elements are sorted. 
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Diagram 2.9.1 Sorting a list by algorithm mergesort. 
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2.10 Parsing 


At this point we are almost able to construct a simple CAS. The only rather 
complicated missing link is an algorithm that is able to recognize legal 
expressions and to decompose them in a form suitable for further processing. 

Two general approaches are used for this operation called parsing: 
bottom—up and top-down (or: recursive descent). Bottom—up methods put pieces 
of the input together until a legal expression is constructed. Top-down 
methods look for a legal expression by first looking for parts of a legal 
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expression, then parts of parts, etc., until the pieces are small enough to 
match the input directly. ; 

In this section an example will be given of a bottom-up parser able to 
recognize legal expressions. A legal expression is constituted by a 
context-free grammar. The context-free grammar defining the set of all 
expressions which should be considered legal to our parser, is given in table 
2.10.1. 


Table 2.10.1 A context-free grammar. 


<expression> ::= <term> | <term> <adding operator> <expression> | 
<adding operator> <term> 
<term> i= <factor> | <term> <factor> | <term> * <factor> | 


_ <number> | <term> * <number> 


<factor> t= (<expression>) | v 
<number> = d | d <number> \ a. # 
<adding operator> ::= + | - f 


a) 


© 


This grammar describes arithmetic expressions such as —7*a(3b+4ab)(a—b). Each 
line is called a production rule. Productions consist of terminal symbols (, 
), +, -, *, non—terminal symbols <expression>, <term>, <factor>, <number>, 
<adding operator> which are internal to the grammar, and metasymbols ::= ("is 
a”), and | (”or”). A v means any character in the range ’A’ .. °Z’, a d means 
any digit; these symbols are also in the category of terminal symbols. 

Below four algorithms are given based on the syntax above: combined they 
are able to check if an expression is legal. In these algorithms the variable 
sym is a global variable and depth (a counter mdicating the recursive depth) 
is initialized to zero. 


Algorithm 2.10.1 Read ~-expression. 


1. [Get character] sym <— next character and increase depth with one. 

2. [Adding operator?] If sym is not an <adding operator> then Read—term. 

3. [Loop] If sym is an <adding operator> then set sym <— next character, call 
algorithm Read—term and after returning, repeat this step. 

4, [End reached?] Decrease depth by one; if depth is zero and the end of the 
input is not reached then give an error. 
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Algorithm 2.10.2 Read—term. 


1. [Number or variable?] If sym is a digit then Read—number otherwise 
Read — factor. 

2. [Test] If sym is not a ’(’, ’* or a variable-identifier then finish. 

3. [Remove ’*’] If sym is a ’*’ then sym < next character. 

4. [Loop] Go to step 1. 


Algorithm 2.10.3 Read - factor. 


1. If sym is a not a ’(’ then go to step 2 otherwise call Read—expression. 
After returning, if sym is a ’)’ then go to step 3 otherwise give an error. 
2. If sym is not a variable identifier then give an error. 


3. sym <— next character. Finish. 
Algorithm 2.10.4 Read-—number. 


1. sym < next character. 
2. If sym is a digit then go to step 1 otherwise finish. 


In these algorithms there is only one recursive call: in Read-—factor the 
algorithm (procedure) Read- expression is called. 
As an illustration of the successive calls to the four algorithms when 


testing an expression, consider the two examples below. 


Example 2.10.1 (a) (b) 
Parse 3a+a(4+c)b Parse 3ab7 
Read — expression Read — expression 
Read —term Read —-term 
Read—number 3 Read-—number 3 
Read - factor a Read — factor a 
+ Read — factor b 
Read — term error.7 
Read - factor a 
Read — factor 


Read - expression 
Read —term 
Read-number 4 
+ 
Read —term 
Read —factor d 


) 
Read - factor b 
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Up to this point we are only able to check if an expression is syntactically 
correct, but we would also like to store a correct expression for further 
manipulation purposes. This can be done by storing the expressions 
unnormalized in an expression tree, which can be done rather fast. In our 
case, we wish to. store a correct expression in a data structure formed by 
expressionptr. Since by using this data structure we can only store 
expressions in expanded form, we have to multiply expressions with each other 
during parsing. After a sub-expression has been read, it will be multiplied 
with the expanded expression which is a result of previous multiplications. 
This will give as a result again an expanded expression. Usually this method 
should not be preferred: if an incorrect expression has been entered, then it 
might take some time before the user is notified of this fact. Consider for 
example the syntactically incorrect expression (a+b)(c+d)(e+f)(g+h)(i+)) 
(k+l). Only after performing five ”expression” multiplications, the syntax 
error + will be recognized. In some cases however this method can be 
acceptable. For example, if it is known that the input expressions do not 
contain many parentheses. Since our goal is to build a simple CAS, we ignore 
the fact that detecting a syntax error might be slow.’ 


2.11 VSCAS - a Very Simple Computer Algebra System 


Using all the information as stated in the previous sections, we are now 
almost capable of building a simple CAS. Only a few algorithms are missing. 
Two algorithms that are able to complete the normalization of an expression 
have not been discussed: an algorithm that adds the equal terms and an 
algorithm that erases the zero terms. And, finally, an algorithm that writes 
the resulting expression on screen. These will not be discussed, since they 
are straightforward. 

By translating all the algorithms discussed in this chapter in Pascal and 
selecting the appropriate data types and declaring the necessary variables, we 
can get an outline of a Very Simple Computer Algebra System in Pascal as 
indicated in diagram 2.11.1, where your-expression is of type expressionptr. 

This CAS is capable of manipulating one-letter characters as variable names 
and unlimited” sized integer numbers. The operatidbns defined on these objects 
are addition, subtraction and multiplication. As has been indicated in the 
previous section, the procedure test—read—and-—store is an implementation of a 


parser combined with multiplication in order to be able to store the 
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expressions in expanded form in a data structure of expressionptr. 

VSCAS is an example indicating that there is no magic behind computer 
algebra. VSCAS per se is not very useful, but it can serve as a basis for 
developing real, i.e. useful, computer algebra programs. Indeed the main 
algorithms of many computer algebra programs discussed in this thesis, ERA for 
example, are derived from VSCAS. In the Appendix the complete Pascal source of 
VSCAS can be found. 


Diagram 2.11.1 Outline of a Very Simple Computer Algebra System. 


begin 
initialize; 
while not end-session do 
begin 
read — test — and — store - expression(your — expression); 
if not error then 
begin 
Sort — expression(your — expression) ; 
add — equal —terms(your — expression); 
erase — zero —terms —in - expression(your — expression) ; 
write —expression(your — expression) ; 
end; 
end; 


end. 


2.12 Conclusion 


The purpose of this chapter was to give an indication of what possibilities 
are available for algebraic formula—manipulation. A short introduction should 
have made the reader aware of the existence of essentially easy to use 
computer algebra systems. However, since many people — even those with much 
programming experience - would not have the faintest notion how to manipulate 
a formula by computer using a higher-level programming language, many sections 
of this chapter have been dedicated to this problem. Algorithms have been 
given which can be implemented using a ”conventional” programming language 
like Pascal or C. Lisp could also be the basis of computer algebra programs. 
In fact, the well-known CAS Reduce is based on Lisp. Lisp-based programs are 
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dependent on Lisp-interpreters. Computer algebra programs using Lisp are 
rather slow compared to those using a higher-level programming language. This 
was one of the reasons why the Waterloo—group used the programming language C 
to create Maple (1983). 

In the rest of this thesis computer programs will be discussed that are 
constructed using the basic algorithms given in this chapter. Especially ERA, 
which is a program that computes the null-space and rank of a matrix 
consisting of algebraic expressions, is a combination of the algorithms given 
in dpis chapter, combined with one more algorithm. 


ry 


2.4 Appendix 
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CHAPTER 3 


COMPUTERIZED EVALUATION OF RANK CONDITIONS: 
IDENTIFICATION OF SIMULTANEOUS EQUATION MODELS 
WITH MEASUREMENT ERROR 


3.1 Introduction 


Studies on the identification situation of structural models often present a 
rank condition, based on the Jacobian matrix of the moment equations, which is 
necessary (under regularity) and sufficient for the local identification of 
the model parameters (cf. Wegge, 1965, Hsiao, 1983, Aigner et al., 1984, 
Bekker and Pollock, 1986). Usually such conditions are rather complicated and 
one is often left with the question how to check the rank of the relevant 
Jacobian matrix in practice. Therefore, most authors also provide order 
conditions or counting rules, such as assignment conditions (cf. Geraci, 
1976), being relatively simple necessary conditions for identification. 

In section 3.2 we will discuss some basic properties of the identification 
problem. In section 3.3 we provide a practical tool for evaluating rank 
conditions for identification in a variety of models. Based on a method 
described in Bekker (1987), a program has been implemented (Exact Rank 
Analyzer: ERA) that uses computerized symbol manipulation to compute the exact 
rank of rather general Jacobian matrices. 

In this chapter — based on a paper by Merckens and Bekker (1988) — this 
approach will be applied to the identification of simultaneous equation models 
with measurement error. For this important class of models another program has 
been implemented - ERASIM - that uses the parameter restrictions as input in 
order to compute the rank of the Jacobian matrix, which itself is also 
generated by the program. The output characterizes the identification 
situation of the individual parameters of the model. 

The formulation of the rank condition for identification in these 
error-shock models used in ERASIM is essentially different from the ones used 
by Geraci (1976) and Hsiao (1976, 1983). It appears that these earlier 
published results are not well-suited for evaluating the identification 
situation in the model. That is to say, when computing the ranks in these 
conditions, problems arise that are not easy to solve. The difficulties are 
due to the fact that the Jacobian matrices are formulated in terms of the 


covariance matrix of the observed endogenous and exogenous variables. Given 
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the covariance matrix of the measurement errors, this covariance matrix may be 
restricted as a result of the model specification, and it is not clear at all 
how to use these non-transparent restrictions when evaluating the rank 
conditions. Consequently, as will be shown, ignoring the restrictions may lead 
to the wrong conclusion. As a result, the earlier published rank conditions 
are useful only to formulate necessary order conditions for identification. . 

In section 3.4 a rank condition will be presented that circumvents these 
difficulties. The exact rank in the condition can be computed by ERA or, 
implicitly, by ERASIM. Section 3.5 gives two examples, and section 3.6 


concludes. 


3.2 Basic concepts of identification 


The classic rank condition for the identification situation in a system of 
simultaneous equations can be formulated as follows. Consider a system of 


simultaneous equations: 
By; = ['x; + €3, j=1,...,t, (3.2.1) 


where B is a non-singular mxm-matrix, and is an mxk-matrix. The vectors x; 
are fixed and the disturbances ¢; are independently, normally distributed with 
zero mean and variance E(e,e; )=2. The parameter matrices B and I are linearly 


restricted according to 
vec'{(B,P)'}R = 1, (3.2.2) 


where both the vector r and the matrix R are known. A well known result (cf. 
Koopmans and Rubin, 1950) says, that a necessary and sufficient condition for 
the global identification of the structural parameters of the i-th equation, 


under a normalization of the diagonal elements of B is given by 

p((B,P)R,) = m-1, . (3.2.3) 
where p(.) denotes the rank, and R; corresponds to the restrictions on the 
i-th equation: the i-th row of (B, I). Model (3.2.1) is identified if all 
equations are identified. 


The result above was found by Koopmans and Rubin using an ad hoc method. We 
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will now consider a more general method to consider the identification 
situation of the parameters of a model. Consider a structure S that describes 
the probability distribution function P(y;9) of a random vector y as a 
function of an [-dimensional parameter vector 6, where 0 € S, and S is an open 
subset of R’. The set of all a priori possible structures is called a model. 
That is to say that the model is defined by the set {P(y;@)| @ e S}. Submodels 
{P(y;9)| @ © T} are defined by sets of structures T, which are subsets of S. 
We assume that y is generated by a known parametric probability function P(.), 
conditional on @. Hence, .a structure is described by a parametric point 6, and 
a model is a set of points T c R'. So the problem of distinguishing between 
structures is reduced to the problem of distinguishing between parameter 


points. 


Definition 3.2.1 

The parameter points 6, € S, and 6, € S, (or the structures S, and S,) are 
said to be observationally equivalent if they generate the same distribution 
for y, ie. if P(y;4,) = P(y;8,) for all y. 


Definition 3.2.2 

The element OF of the parameter vector 6° is said to be locally 
identifiable (in S) if there exists an open neighborhood of 6° containing no 
point (9 € S), with & #4 a, which is observationally equivalent to 6°. 


If all elements of 9° are locally identifiable then we say that 6° itself is 
locally identifiable. If the open neighborhood in Definition 3.2.2 can be 
taken to be R,, then the identifiability is said to be global. 


Definition 3.2.3 

Let M(@) be a continuous matrix function of 8 € R', and let 6° < R'. Then 
6° is said to be a regular point of M(@), if the rank of M(@) is constant for 
points in R! in an open neighborhood of 6°. 


We now make an assumption that is justified in most practical models. Let 
o(@) be an analytical function a: R' > R", so that there exists a one-to-one 
relation between the elements of {P(y;0)| 6 « R'} and {o(@)| Oe R'}. For 
example, if y is normally distributed, then all ‘information necessary for the 
evaluation of the identification situation of the parameters is contained in 
the first and second-order moments. In this and subsequent chapters we make 
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this assumption of normality, so that the first and second-order moments can 
be considered as a function o(8), as described above. Consequently two 
structures are observationally equivalent if they produce identical first and 
second-order moments. The model specification may restrict the set of 
structural parameters: we assume @ to satisfy R(@) = 0, where R(.) is an 
r-vector of analytical functions. Following Fisher (1966, Theorem 5.9.2), we 


can give the following theorem. 


Theorem 3.2.1 
Let J(@) be the (n+r)xm Jacobian matrix formed by taking partial 


derivatives of | of) with respect to 6, 
R(@) . 
80(@)/00" 
J(0) = ‘ - (3.2.4) 
OR(8)/30" , 


If 0° is a regular point of J(@), then a necessary and sufficient condition 
for 6° to be locally identifiable is that J(@) has rank m at 6°. 


Theorem 3.2.1 forms the basis of all rank conditions, derived in this thesis. 
In order to characterize the identification situation of a model, the exact 
rank of an ‘analytical Jacobian matrix J(@) (3.2.4) should be computed for 
regular points. Furthermore, if J(8°) is net of full column rank, it may still 
be the case that single elements of 6° are identified. If the i-th element of 
0° is identified, then under suitable regularity conditions, the i-th column 
of J(8°) must be linearly independent of the remaining columns of J(6°) (cf. 
Bekker and Pollock, 1986), so that vectors in the null-space of J(0°) must 
have their i-th element equal to zero. 

A description of a computer program for the computation of the exact rank 
of J(@) for regular points, and for determining the linearly independent 
columns of J(@) in case of underidentification, will be given in the next 


section. 


é 


3.3 Computerized evaluation 


As the elements of the matrix function J(@) in (3.2.4) have, in general, no 


fixed numerical values, standard computer programs are not well-suited for 
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evaluating the rank of J(@). In fact, the rank of J(@) is a function which 
varies across the parameter space. However, as was noted by Bekker (1987), the 
rank of J(@) is constant for all regular points, that is for almost all 
parameter values. As the assumption of regularity has already been made in the 
rank conditions for identification, a computer program should just compute the 
rank of J(@) for a regular point. For this problem computer algebra is 
particularly well-suited. That is to say, using symbol manipulation, it is 
possible to compute the exact rank of J(#) for regular points. Below we 
briefly discuss the computer program ERA, the Exact Rank Analyzer, which 
computes the rank of matrices that may have rather general parameterizations. 
In fact it computes a parameterization of the null-space of the matrix, thus 
giving information about which columns of the matrix are linearly independent 
of the remaining columns, which is, as we saw, information relevant to the 
identification of the separate parameters. The program ERA will be used 
throughout this thesis, as a tool for evaluating the identification situation 


in a variety of models. 
3.3.1 The Exact Rank Analyzer (ERA) 


Let A(@) be an mxn-—matrix, where the elements of A(@) are analytical functions 
of a parameter vector 9 of order k. Bekker (1987) has given a scheme for the 
computation of a matrix N(A(@)) which forms a basis of the null-space of A(@) 
for almost all parameter values. Consequently, as was proven by Bekker (1987), 
for any regular point A(@) the rank of A(@) equals n minus the number of 
columns of N(A(@)). Furthermore, zero rows in MN(A(@)) indicate linearly 
independent columns of A(@). The matrix N(A(@)) is defined recursively as 


follows. 


(a) If m=1, so A(@) is a row vector, A=a’ say, then three cases are 
distinguished for N,=N(a’): 
(i) if @’=0 id. then N,=/, 
(ii) if n=1 and a’#0, then N,=0 
(iti) if n>1 and a@’#0, let a’e; be the first non-zero element of a’, 
where e; is the i-th unit vector, then 


N,=(e,a' -(@ ep) Tn) (€1y--y€¢-1)€:419-9En) - 


(b) If m>1, let A=(a,,...,a,,)’ then 
Ny = N{(Gq,--501)' } = N{(Gy,.-541-3)’} N[ ay’ N{(G1,.-,0;-1)'}), 
for l=2,..,m. 
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For a more complete description of the matrix function N(A(9)) and _ its 
properties see Bekker (1987). This scheme can easily be computerized. The 
program ERA performs this recursive ir on a matrix A(9)y, where its 


elements are assumed to be analytical expressions: of the form: 
8 k b5; i : 
>» C3 Il 93; 3 # (3.3.1) 


where the c; and Pigs t=1,..,5, and j=l,..,k, are integers while the bj; are 


non-negative. , . 
For example, consider the 6x5—matrix A, parameterized with 10 parameters 


(a, b, c, d, e, f, g, h, i, 9): ; Sea 2 
: 
j d 00 j 
i e 0 Ore 
a+d 0 be b 
A= h b+e f g h ” 
e a+b f ge a 
d*>+hi 0 0 0 d°+hi : 


For N,,...,V_ we find: 


| @ 0 0 j 0 01 
-j 0 0 0 0 0 0 
M=| 0-j 0 O|, N= j(ej-di) |-1 0 0 |," 
0 O0-j O "0-1 0 
0 0 0O-j 0 0-1 
. ~“ 
0 b bg-cf 
0 0 0 
Ns = j(ej-di)’ | ¢ -at+b-d |, Ng=N5=Ne= bj*(cj-di)* | —ag+bg—dg 
-b 0 af ~bf+df 
0 -b 5 -| -bo+cf 


. 
x 


Consequently, for all regular points the rank of A equals 5-1=4. The zero on 
the second position of Ng indicates that column 2 of matrix A is independent 
of all,the other columns. Actually ERA gives as output Ng without the factor 
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bj*(ej-di)*. A part of this blow-up’ factor is systematic: at each step it is 
the first non-zero element of a;’N(q,,..,¢,,)', where a; is the last row of 
(@1,-.,@,4)’for which  a;'N(q,,..,q_,)'#0. ERA makes these important 
systematic simplifications at each step. Furthermore, it checks, at each step, 
to see whether further simplifications can be achieved. These checks do not 
exclude the possibility that the output — the null-space parameterization — 


can be simplified even further. 


3.4 Contemporaneous error-shock models 


In the Handbook of Econometrics (1983, Vol.1, pp.257~265) Hsiao considers, 
among many other ‘problems, the identification of the following simultaneous 


equation model with errors in variables. Let an economic theory specify 


Bn t+ Ta = %, 
v= t+ &, (3.4.1) 
X= Oe + &, 


where 7, is an unobservable m-vector of endogenous variables, £, is an 
unobservable, stochastic, n—vector of exogenous variables, and (, is a 
m-vector of disturbances; B and [ are mxm and mxn constant matrices, 
respectively, and B is assumed to be non-singular. The vector of disturbances 
¢ has zero mean and is uncorrelated with €,; its covariance matrix is denoted 
by Z. 

Instead of the latent variables m, and £, we observe y, and x, 
respectively. The error vectors €, and 6, have zero means and covariance 
matrices @, and 5, respectively, which may be singular. The vector 6, is 
assumed to be uncorrelated with €,, 7, and €,; the vector €, is assumed to be 
uncorrelated with &, (Hsiao uses this assumption, but he does not mention it) 
and 7, (this last assumption is not really necessary for the analysis given 
below). All variables are assumed to be identically and independently 
distributed over t. The covariance matrix of €, denoted by Cg is 
unrestricted. As we shall use only second-order moments to solve the 
identification problem, we could as well make the assumption that all 
variables are normally distributed, which is a conservative assumption when 
dealing with problems of identification (cf. Aigner et al., 1984, Bekker, 
1986). 
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Rewriting model (3.4.1) in terms of observables, we have 


By, + Tx, = ( + Bez + T&, (3.4.2) 


. 


and thus all (second-order) observational information is contained in, 


. ; , Re at 
Cy = E(x, x) = EE, & ) + 5 = Ce + Os, + (3.4.3) 
Cy = Ely % ) = BTC, - O5)I"B™ + BE BT + &, ; 
= BT Cer'B +E, ~ (3.4.4) 
Cy = lye xt) = -B'T(C, = 5). (3.4.5) 
.# 


When © and &, are unrestricted, these matrices cannot be identifted from these 
equations. However, this does not affect the identification of B, I arid. B;..- 
Therefore, ae C. + Be, will be treated as an amalgamated disturbance term 
with covariance matrix 5”. If B, T and @s are identified from these equations, 
so will be 5". 6 

For the identification of B, [ and %s; we need only consider: equation 
(3.4.5) together with the restrictions on B, [ and Ss (which we assume to be 
linear for presentational convenience; the analysis can be easily generalized 
to include non-linear restrictions as well, cf. Bekker, 1987): 


Rwe,r) heel =d, (3.4.6) 
vec(I") 
Rg, vec(®s) = 6. ' (3.4.7) 


The constant matrices Rig r) and Ro; are of order T(B,r)yX M(m+n), and Ty 5% n?, 
respectively; the constant vectors d and e-are of order rg y)x 1, and rg 5% 4; 
respectively. The number of restrictions on B and I is given by tgp), and 
Tbs denotes the number of restrictions on ,, including the %n(n-1) symmetry 
restrictions on ®;. Equation (3.4.5) can be written, after pre—multiplication 


and transposing, in stacked form as 
vec{C,,B'+ (C,-®5)I"} = 


vad| aa) 


[Un @ Cyx)) (Im @ (Cx—-5)) | laa 
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Thus the matrix of derivatives of (3.4.6), (3.4.7) and (3.4.8) with respect to 
the elements of B, F and 5, that is with respect to 0’ = (vec’(B’), vec’(I"') 
vec'(®s)), is given by 


vec{C},B' + (Cz-%5)I"} 


ies 


} / 00" = 


Rg ,vec(®s ) 


(Im ® —(Cx-O5)P'B™) (Im ® (Cy-O5)) (“I @ In) 


iw 0 
Re,ry | ™oo™ Rw,r) 0 . (3.4.9) 
0 In ® In 


0 0 Re, 


Consequently, restating Theorem 3.2.1 we find the following result. 


Theorem 3.4.1 
If 6° is a regular point of J(@), then a necessary and sufficient condition 
for (B, IT, 5, ry to be locally identified is that J(6°) has full column 


rank. 


With respect to an analogous Jacobian matrix, which will be discussed below, 
Hsiao (1983, p.261) made the statement that, using an instrumental variable 
argument, it can be shown that if the structure S is locally identified, it is 
also globally identified. However, we have trouble finding the argument and, 
in fact, we are doubtful about the statement. The reason is that the equations 
in (3.4.8) are not linear in the parameters of B, [ and 9%;, so that the 
possibility of two locally unique solutions to the equations in (3.4.6), 
(3.4.7) and (3.4.8) cannot be excluded. 

Before discussing Hsiao’s Jacobian matrix we first want to simplify the 
rank condition in Theorem 3.4.1. Making the common assumption that the matrix 
Cg = C,-G, is non-singular, we premultiply J(@) by 
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—(Im @ Ce )* 0 0 


vs a Trp ry & (3.410) * 
0 0 I : 
"05 j 
and postmultiply by 
(Im ® B’) 0 0 . # 
W= |(In @I") (Im @ In) (@ Ip )I- (3.4.11) 
0 - 0 (In ® Ce) ee 


As both matrices V and W are non-singular, the resulting matrix yaw will be of 
full column rank if and only if J is of full column rank. Furthermore, the 


matrices J and VJW will share regular points. We find 


0 - (Un ® In) 0 
, In ® B’ 0 0 
= |Rery)| ™ Rip,r ) Rg,r) . 
VW = Ler eae reln .* (3.4.12) 
.lo¥ 
O° 0 Ros (In ® Ce) 
Defining 
In ® B _| 0 
% Ren| aie | Re, o| | 
J(9) = In @ T rel, ; (3.4.13) 
0 Ro (In ® C¢) ~ 


we find that rank(J(@)) = rank(J(8))+mn. Consequently, J(@) is of full column 
rank if and only if J(9) is, and J(@) and J(@) share regular points. 


Corollary 3.4.1 
If 6° isa regular point of J(8), then a necessary and sufficient condition 
for (B, I, 9s, 5" to be locally identified is that J(0°) has full column 


rank. 





The rank condition in this result is transparent in the sense that if ®;=0 — 


no measurement errors in x — then Ro, is equal to J,2 and consequently 


¥ 
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non-singular, so J(8) is of full column rank if and only if Rypr) (Um ® B), 
(Im @ [))' is of full column rank. The latter condition corresponds to the 
classical condition for identification in a simultaneous equation model 
without errors (cf. Schmidt, 1976). 

It follows from Corollary 3.4.1 that a necessary condition for 
identification is that J(8) has at least as many rows as it has columns. So an 


order condition for identification is given by 
T(B,r) + "$5 > m24n2. (3.4.14) 


As rT ; includes we symmetry restrictions on $5, we may also write ry 5 = r% * 
yon(n—1), where rg . corresponds to the number of restrictions on the diagonal 


and subdiagonal elements of $5. Consequently, 
65 > m2 Tp, r) + Yin(n+l) (3.4.15) 


gives a bound on the number of restrictions on the variances and covariances 
of the errors, which should be satisfied for the model to be identified. Of 
course such an order condition is not sufficient for identification. 

In order to characterize the identification situation of a model, the exact 
rank of 70°) should be computed. Furthermore, if 76°) is not of full column 
rank, it may still be the case that single elements of 6° are identified. If 
the i-th element of 6° is identified, then the i-th column of J(0°) must be 
linearly independent of the remaining columns of J(0°) (cf. Bekker and 
Pollock, 1986), so that vectors in the null-space of J(9°) must have their 
i-th element equal to zero. A description of a computer program for the 
computation of the exact rank of J(@) for regular points, and for determining 
the linearly independent columns of J(@) in case of underidentification, will 
be given in the next section. Section 3.5 gives examples. 

First we shall briefly discuss the Jacobian matrix derived by Hsiao (1983). 
He took derivatives of (3.4.6) and (3.4.8) with respect to B, [ and only the 
free elements in 5. For the characterization of the identification situation 


in several examples he evaluated the rank of the following Jacobian matrix: 
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(Im ® Cyx), (Im @ (C,-@5))  U 


. aoe. 0 : 
Rar! "°°" Ru,ry 0 
0 In ® In 


where, in the present notation, U = (-I @ I,)L, and Rg 5 = 0, while (Ry yh) is 
non-singular. The main difference with J(@) is not due to the use of the 


* 


J (3.4.16) 


matrix U; instead it is due to the fact that the matrix C,, has not been 
written as -B“I(C,-®5). This is an important difference since the matrix B'P 
may be restricted due to the restrictions on B and I’. Consequently, the 
elements of C,, may be related to the elements of C,-;. These restrictions 
* may exist whether or not the model is identified (cf. Fisher, 1966, Bekker and 
Dijkstra, 1990). However, the relations between Cy, and C,-%s are far from 
being transparent in the Jacobian matrix J”. As a result the evaluation of the 
rank of J” (or similar matrices presented by Geraci, 1976, or Hsiao, 1976) may 
lead to confusing conclusions. For example, in section 3.5a model will be 
specified for which J” is of full column rank for almost all values of C,,, 
thereby excluding only sets to the Lebesgue measure zero, while J(@) is of a 
deficient rank for all regular points so that the model is, in fact, 
underidentified, contrary to what an evaluation of Cd might suggest. 

As we have shown in this section, the evaluation of the Jacobian matrix in 
(3.4.13) gives a correct characterization of the identification situation in 
simultaneous equation models with errors in the variables,” like model (3.4.1). 
Since the construction and evaluation of such a Jacobian matrix by hand would 
be rather tiresome, we have created a computer program, ERASIM, which uses as 


input the restrictions on the matrices B, [ and @s to generate the matrix 


J(9), as in (3.4.13). This Jacobian matrix will then serve as input for 


evaluation by ERA. If J(@) is of full column rank for regular points, the 
parameters in B, F and 5 are docally identified. If not, ERA gives back 
N(J(8)), and the second part of ERASIM derives a parameterization of N(J(@)), 
where J(#) is the Jacobian matrix in (3.4.9). This output, N(J(@)), can be 
separated into parts corresponding to the matrices B, IF and 95; it is 
indicative of the identification situation of the separate parameters in the 


model. 
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3.5. Examples 


In this section we analyze the identification situation of the following 


model: 


Model 3.5.1 


Ya= “Yuka ro 

Y= ~Ya18ea-7 228 t2 eri 
¥e3=—-PaVer -Ya2€ 12 esi 
Vea=—ParVer ~Yastis eas 
X=Ee1 s o 0 0 
X=Et2t O12 5 @5 = E(5~5;) =|0 Ago Aza}. 
Xy=Er13 + 513 , 0 Azz A33 


In this model t. = (; + Be, is an amalgamated disturbance term. That is to 
say, y, may or may not have measurement errors, but =", which is the 
covariance matrix of a is assumed to be unrestricted. Only x, and x3 are 
subject to measurement errors, and the covariance matrix of 6. and 5:3 is 
unrestricted. Clearly the first equation is identified, that is Yun is 
identified. The question is whether the remaining parameters are identified. 

It will be shown that the construction of J(8) is a rather tedious job, 
however this job is taken care of by the program ERASIM. This program also 
shows that 7, is the only parameter that is identified. The model as a whole 
is not identified. Furthermore, Hsiao’s J” will be constructed and it appears 
that this Jacobian matrix J” is of full column rank for almost all values of 
C 


yx 
underidentified. Apparently, given the covariance matrix of the measurement 


This might lead to the wrong conclusion that the model is not 


errors, the elements of the covariance matrix of y and x are restricted by the 
model specification. 
In addition to Model 3.5.1 we also consider another, less restricted, 


model: 
Model 3.5.2 


The same model as Model 3.5.1 except for the restriction ,;=0. This 


restriction is relaxed and the first equation becomes: 


* 
Ya=—-Y¥uSa —Vi36e3+¢ e1- 
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It appears that this model, with a less restricted set of parameters, is not 
underidentified. For all regular points of the Jacobian matrix J(6), that is 
for almost all parameter values, Model 3.5.2 is identified. The extra 
restriction 7;3=0, which holds only for a set of Lebesgue measure zero, leads 
to Model 3.5.1 and it renders the model underidentified. 

For Model 3.5.1 we have the following parameter matrices: 


1000 Yu 9 O 0 0 0 
0 140 

B = Of, ra] M29 |, be =] 0 Age Avs 
3,9 1 0 0 32 0 O Aes des 
Ba 0 0 1 0 0 ¥43 


If we consider the necessary order condition (3.4.15) for identification, we 
Fy 
find that , = 3, m? = 16, gr) = 21, and ¥n(nt1) = 6. So 


* 


Tos = 32m - rer) t+ ¥n(nt1) = 1 


is satisfied; there appear to be two ’overidentifying’ restrictions, and the 
model may well be identified. The input of the ERASIM program consists of the 


information: 8,=1, 62=0, 83:=Ps, Ba=Pa, Pi2=0, ~~) Yu=Yuy Ya=721, 
Ya1=0, 5 A=O, Agr=0, Agi=O, Aze=Azz, Azg=Azay Aga=Azg- Using this 
information it constructs the Jacobian matrix J(@). As the elements of B and I 


are not linearly related, the matrix Rig rp) can be separated in two parts: 


Rg, 0 
Rar) =| 7 
0 Rr, 
where Rp'vec(B’) = dg , and Rp'vec(I") = dp , say. The Jacobian matrix J(9) 
can thus be written as 
Rz'(Im @ B’) 0 
J(0) =| Rr'Um @ I") Rp'(F ® Ip) | (see 3.4.13). 


0 Rg,(In ® Ce) 


For Model 3.5.1 the submatrices of J(@) take the following form. 
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1 0 B3;Ba1 
o 1 0 6 0) 
00 1 °0 
00 01 
1 0 B31841 
0100 
; i oo 10 
Rz'(In ® B’) = 0001 ’ 
0100 
4 0010 
0001 
16) 0100 
0010 
0001 
0 Y22 Y32 0 0 
0 0 0 43 
0 0 0 Y43 
Rr'(UIn @ I’) = Yu Ya 9 ’ 
0 0 0 43 
O Y11 Ya 9 0 
8 Y22 Ys2 0 
071° 000000 
0 071, 9 6 0 0 0 O 
0 OY, 0 0 Y22 0 0 0 
Rr'(l @ In) = 0 0 0320 0 0 0 of, 
0 0 0 0 0 Y¥32 © 8 0 
0 0 0 0 0 043 0 0 
0 0 0 00 0 0743 0 


Cig Cop Cop 09 0 08 0 80 OO 
Cig Cop C93 0 20 0 0 0 0 

Ro sUn ® Cg) = ’ 
0 0 0 Cy C42 43 0 «208 O 


0 0 0 0 0 DO Cy CyQ C43 


Oo 0 0 43 C23 €33 -Cy2 -C22 -C23 


where Cg = C, — 5 has typical elements ¢;;- It should be noted that apart 

from the symmetry, the elements of Cg are unrestricted. Consequently, the 

diagonal and subdiagonal elements of C, can be treated as free parameters. 
Obviously, the construction of F(8) is rather tedious, however this 


construction has been taken care of by the program ERASIM, which also computes 
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a parameterization N(J(@)) of the null~space of J(@). In case of a non-trivial 
null-space it will subsequently derive a parameterization of N(J(@)), the 
null-space of J(@). This latter null-space is informative with respect to 


which elements of 9 are underidentified; it can be found by the operation 


(In ® B’) 0 
N(J(8)) = |Um @T’) (P@® In) | N(I(8)). 
0 (In ® Ce) 


For Model 3.5.1 ERASIM gives a one-dimensional null-space of J(@), so that 
’J(@) has a deficient column rank for all regular points, and Model 3.5.1 is 
consequently underidentified by Theorem 3.4.1. The null—vector, n(@) say, 
satisfying the equations in (3.4.6), (3.4.7) and (3.4.8), which is found by 
ERASIM, can be presented as follows. 
vec( Ag) 
n(@) = | vec(Ap) |, 
vec(Ay 3) 


where the elements of vec(Ag), vec(Ay) and vec(Ay ga correspond to the columns 
of J(@), with the derivatives taken with respect to vec(B’), vec(I”’) and 


vec(®s), respectively. We find 


0 00 0 0 0 0 
0 00 0 Cy2¥22 —C117: 0 
Ags An i= 12%22 —€11Y22 
—C12¥32 09 0 O |’ a ie —€11 732 0 ; 
C3743 0 0 0 0 —€117 43 
e 0 0 0 
2 
Ag = oan O  C42-€4 1 C22 Cy2€ 13—€11 C23 


2 
O €y2€y3-C€41€23  €13—€1 1 C33 


Ag shows the underidentification of 83, and 84;;, Ar shows the identification 
of Yu, and the underidentification of all remaining free elements of [; Ay P 
shows the underidentification of Az2, Agog and Ag3. 


When we consider Hsiao’s Jacobian matrix 
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(Im ® Cy) Um @ Cg) U 


J = Ry’ 0 0 | (see 3.4.16), 


0 Rr’ 0 


where U = (-I @ I,,)L, and Ry L = 0, while (Rj_,L) is non-singular, we find 
5 o5 § 


$11512513 514 
$21522523 5 24 
531532533 S34 Oo 
$11812513514 

$21522523524 

In ®Cxy= 531532533534 


$11812513514 
$21522523524 
531532533534 


$ 11512513514 
O S$ 21522523524 
531532533534 


where Cy, has typical elements s,;. The matrices Ry " and L are given by: 


oocsc 2 


-1 


oH 
eococcow 
eoocore 
coorono 
coreoe 
cocooce 
or oooo 
ke ocaco$o 
eooaooe 
coocoooo 

- 

S 


so, consequently, the matrices I, @ Cg and U = 


€31€12€13 
C1222 6 23 0 
C13C23 C33 
€41€42€13 
Cy2€92C23 
In ® Ce = eae ; 
11€12013 
€12€22C23 
€13C23C33 
0 © 1101213 
€12€22C23 
€13€23C33 


Finally, Ry’ and Ry’ are given by 
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ooo00100000 
oo0o00001000 Ff, 
qooo00000001 


(- @ I,,)L are given by 


0 0 0 

0 0 0 

0 o oO 

0 o 0 
Yo2 9 0 
0 Y22 0 

,U=- 22 

0 0 0 
Y32 9 0 
9 Y32 6 

0 v o 

9 Y43 0 

{ 8 0 Y43 





ge 
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_ 


1000 
0100 O 
0010 
0001 
1000 010 0 
0100 ; ceo 
As! 0010 , 
Rp'= 0001 » Rr'= 
0100 
0010 O aa 
0001 ; fees 
0 0100 
0010 
0001 


Evaluation of the rank of z; by hand or, preferably, by the Exact Rank 
Analyzer, ERA, shows that J” is of full column rank for almost all values of 
Cy, This may lead to the wrong conclusion that Model 3.5.1 is identified. The 
reason for the discrepancy in the ranks of J(@) and I is that the matrix C,, 
is, apparently, restricted by the elements of C,—5. These restrictions on Cyx 
are far from being transparent and it is not clear at all how these | 
restrictions should be taken into account when evaluating the rank of J”. The 
alternative Jacobian matrix J(@) does take the restrictions on Cy into 
account; in fact the matrix Cy, is not present in J(@). Therefore evaluation 
of J(8), using ERASIM, leads to the correct, unambiguous solution to the 
identification problem in simultaneous equation models with errors in the 
variables. a 

For Model 3.5.2 the matrices B and % are subject to the same restrictions © 
as in Model 3.5.1; only I is restricted differently: 


1000 Yu 9 13 00 0 
0100 0 ‘ 
B= P= | Ta Yar" |, By = | 0 Ase Aas 
B31 0 1 0 . 0 Y32 0 0 roa Aas 
01 00 43 


Bai 0 


When constructing J(@) we find that the submatrices Rg’ (Im @ B) and R(I, ® Cg) 
are the same as ‘in Model 3.5.1. The matrix Rr'(I,, @ I) becomes 
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0 Y22 Y32 9 O 
Yizs 9 0 Yaz 
Rr'(Im @ I") = ya hie 
Oo Yu Yai 0 0 
0 ¥22 ¥32 9 


and R-(fel,) is given by 


Yu9 0 © O oO 
Yai 9 +9 Y22 0 


0 

0 07320 0 Oo 
0 0 0 0 Y32 0 
0 
0 


ooc oc 


R(T @ I) = 


0 0 0 0 Y43 0 
000 0 0 Y4g 


oooe os 
oo coco coe 


Evaluation of the identification situation of Model 3.5.2 by using ERASIM 
shows that J(@), and J(@), are of full column rank for all regular points. 
Consequently, Model 3.5.2 is identified. It may be useful to note that 7(8), 
and J(@), are not of full column rank for parameter points for which y,; = 0, 
corresponding to the extra restriction in Model 3.5.1. However, these 
parameter points are no regular points of J(8), or J(@), and the statements of 
Theorem 3.4.1 or Corollary 3.4.1 refer only to regular points. From the 
analysis of Model 3.5.1 we know that if y,;; = 0 then, locally, observationally 
equivalent points do exist; the analysis of Model 3.5.2 says that if y,, #0 
then there are no observationally equivalent points locally. 


3.6 Conclusion 


It has been shown that an evaluation of the identification situation in 
error~shock models using the rank condition presented by Hsiao (1976,1983) is 
hazardous. The same holds, although is has not been shown in this chapter, for 
the rank condition presented by Geraci (1976). As an alternative we have not 
only presented a correct rank condition, but also two computer programs, ERA 
and ERASIM, which take care both of the tedious job of constructing a suitable 
Jacobian matrix and of computing its rank for regular points. It is a 
practical tool for analyzing the problems related to the identification of the 
separate parameters of the model. 

The program ERA is suited for evaluating rather general rank conditions for 


identification. These include the rank conditions for local identification in 
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simultaneous equation models with covariance restrictions on the disturbances 
(cf. Bekker and Pollock, 1986, Hausman et al., 1987) and similar conditions 
for factor models (cf. Bekker, 1987). In chapter 4 the identification problem 
in general restricted factor analysis models will be discussed, where we again 
use ERA, as well as in chapter 6 which treats the: identification problem in‘ 
general LISREL models (cf. Jéreskog and Sérbom, 1984). 
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CHAPTER 4 


COMPUTERIZED EVALUATION OF RANK CONDITIONS: 
IDENTIFICATION OF RESTRICTED FACTOR ANALYSIS MODELS 


4.1 Introduction 


In this chapter a rank condition for the local identification in restricted 
factor analysis models will be presented. This rank condition is similar to 
the rank condition as presented by Bekker (1989): it gives a rank condition 
for the local identification of the restricted factor analysis model as a 
whole. A practical computer program will be presented that uses the parameter 
matrices as input in order to construct the Jacobian matrix, which can be 
evaluated by ERA in order to characterize the identification situation of the 


individual parameters. 


4.2 Restricted Factor Analysis Models 


The general Factor Analysis (FA) model specifies the set of linear 
relationships 


y= AE +e, (4.2.1) 


where y is an observable p-—vector of endogenous variables, ¢ is an 
unobservable, stochastic k-vector of latent exogenous variables (also called 
common factors) and € is a p—vector of disturbances (specific factors). It is 
assumed that the disturbances € are normally distributed and independently of 
g. A is a pxk-matrix of coefficients (factor loadings). The covariance 


structure is thus given by 
= AGA + ¥, (4.2.2) 


with 2, ® and W the covariance matrices of y, € and e respectively. 

In restricted (or confirmatory) FA the matrices 6, A and W are subject to 
linear restrictions and W is assumed to be diagonal. Here model (4.2.2) will 
be considered where $, A and W are subject to arbitrary restrictions. 


Under the assumption of normality, the moment equations of the FA model 
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(4.2.2) contain all information relevant to’ the identification of the 
unrestricted (free) parameters in A, ® and ¥. The vector of these unrestricted 
parameters will be denoted by 8, so the elements in A, ® and W are functions 
of this parameter vector. It is assumed that 6 has a one-to-one relation with - 
the elements of A, ® and ¥. 

The matrix of derivatives of (4.2.2) with respect to @ is given by 


° 


8 vec { A(O) B(8) A'(8) + (A) } 
J(@) = - = 





a 6" ~ es 
= (Ip2 +.Kpp)(Ip @ A) L + (A@A) F + P, (4.2.3) 


with L = 8 vec(A’) / 00', F = @ vec(%) / 060’ and P = @ vec(¥) / 06’, and Ky» 


is the p’xp?- commutation matrix. A commutation matrix Ky,q is defined to be a 


pqxpq-matrix consisting of gp blocks of dimension pxq, where element (j,#) of 
block (i,j) has value one and the rest of the elements of this block have 


value zero. Since @ has a one-one relation with the elements of A, ® and ¥, 


: L 
the matrix | F | has full column rank. Restating Theorem 3.2.1 we get 
P 
Theorem 4.2.1 


If 4 is a regular point of J(@), then a necessary and sufficient condition 
for 99, and thus for the structure of 5 = (A, , W)g, to be locally identified 
is that J(@,) has full column rank. 
Now assume that ¥ is linearly restricted. Also assume that and (A, ®) are 
functions of separate parts of. the parameter vector 9’ =(6:,03), A=A(6,), 
b=6(6,), Y=(0,), say. Then the Jacobian : 


Ry Avec {A(8,) B(A,) A’(8,)} 


J(0,) = (4.2.4) 





30; 


can be considered as well for determining the identification situation of the 
parameters in (4.2.2). The matrix R contains fixed elements and corresponds 


to the linear restrictions on ¥: 


R, vec (W) = Fi Bs (4.2.5) 
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Using some basic matrix calculus it follows that 
J(6,) = R,, {tye + Kpy)(Ip @ AD) L + (A@A) Fh, (4.2.6) 


with L = 8 vec(A’) / 06; and F = 4 vec(®) / 80}. 


Corollary 4.2.1 

If W is linearly restricted according to (4.2.5) and 0, is a regular point 
of J(@,), then a necessary and sufficient condition for the structure S = (A, 
@, W)y to be locally identified is that J() has full column rank. 


Since $(8,) is a symmetric matrix, A(0,) £(6,) A’(9,) will also be symmetric 
and a row in Ry, representing a symmetry restriction on ¥(6,) will result in a 
zero row in the Jacobian matrix. Therefore Ri need not contain the symmetry 


restrictions on W. 


It follows from Corollary 4.2.1 that a necessary condition for identification 
is that J(8) has at least as many rows as it has columns. So an order 
condition for identification is given by 


Ny = 76,5 (4.2.7) 


where nm, corresponds to the number of restrictions on the diagonal and 
subdiagonal elements of ¥ and 79, corresponds to the number of parameters in 
9,: the number of parameters in A and &. 

In order to fully characterize the identification situation of an FA model 
the exact rank of J (9p) should be computed. Furthermore, if J(o) is not of 
full column rank, it may still be the case that separate elements of 6 are 
identified. A necessary and sufficient condition for the i-th element of 8 to 
be identified (for almost all values of 0.) is that the i-th column of J(9y) 
is linearly independent of the remaining columns of J(,), so that vectors in 
the null-space of J (99) must have their i-th element equal to zero. 

A description of a computer program that can construct the Jacobian matrix 
as given in (4.2.6) and analyze its rank and null~space, will be given in the 


next section. 
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4.3 Computerized evaluation 

The Jacobian F(9,) has to be constructed before it can be evaluated by a 

program like ERA (see chapter 3). As the elements of the Jacobian matrix are 

functions of the parameter vector 9, and not just fixed numerical values, the 

Jacobian matrix can only be computed automatically by a computer program that 

manipulates symbols and formulas. Such a program is given by ERAFAC, which 

uses as input the parameter matrices A(@,), $(6,), and ¥(6,) and returns as 

output the Jacobian matrix T(9o) as in (4.2.6). In this section a short 

discussion on the construction of ERAFAC will be given. . 

Formula (4.2.6) could be programmed at face value. However, as is often the 
case with programming formulas and algorithms, a more efficient and simpler 
method can be used. If formula (4.2.6) is used to construct J(6,), then the 
following matrices should be built: Ry, £ and F (which are derivatives of A 
and @, respectively), [@A, A@A and the commutation matrix K,. Then by 
multiplying and adding conform (4.2.6), the Jacobian could be computed. Many 
multiplications and much unnecessary bookkeeping are involved . It is simpler 
to consider (4.2.4) as basis of the computer algebra program and.proceed as 
follows: 7 . 

1. Read the parameter matrices A, @ and W, containing the restrictions of 
model (4.2.2). . , 
Construct R, and AGA’. 

Determine the free parameters 6, in (A,®). 
Take the derivative of AGA’ with respect to 6. 
Premultiply the result of step 4 with R, 


oA Pp wD 


Write the computed result. 


Additionally, a more efficient storage method for Ri can be used: Ry is a 
sparse n,xp? matrix, having just a few elements in each row filled. So instead 
” of storing this matrix in an array it can be stored in a data structure, which 
keeps only information on the non-zero elements of Ri. The example below will 


show the principle of this method. 
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Example 4.3.1 


Wi 
Assume that p=3 and¥Y= |]. . , then 
4p +2 
oder omeud ¢ 
« dys 
2 ee (see 4.2.5), (4.3.1) 

4 ala ~1 

and r, = 0. Instead of storing Ry in an array, it will be stored in a vector 


where each row specifies the same information as the same row in Ry: each row 
contains a linked list (see chapter 2) with the following information: 
(element number, value, next), where elements having zero value are not 
included. So in this example, Ry, could be represented schematically as in 
diagram 4.3.1. 


Diagram 4.3.1 


2 | 


Ri, = 
v 
EES 
EXESE IESE 


The fourth row of diagram 4.3.1, for example, can be read as: the first 
element of the fourth row of Ry, equals 4, the fifth element equals 1, and the 
ninth element equals -1 and the rest of the elements are zero. A data 
structure like the one above is also more efficient in comparison with the 
usual matrix—multiplication method. For example, consider the 
premultiplication of a 9x9-matrix X by Ry. For each element in the resulting 


matrix nine multiplications are needed if a standard matrix—multiplication 
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method is used. A multiplication method that uses the data structure above,’ 
will need less multiplications, because all elements having zero values are 
eliminated in this representation of Ry. E.g. the first element (and also the 
other elements) of the fourth row of the resulting matrix will need three 
multiplications. If Ry, is a sparse matrix indeed, then the memory required for 
storing the (information on) elements of Ry, using a data. structure as in 
diagram 4.3.1, will also be less than the usual array. (matrix) -storage. 


Premultiplying X with R,, gives a matrix of the form: 


Xo, 

x3, 

R'X = ‘ , 
-| Xe. 


4X +X5,-Xoq. 


where x; means the i-th row of X. 

ERAFAC uses the methods described above to construct the-Jacobian matrix 
necessary for evaluating the identification situation in model (4.2.6). The 
resulting Jacobian matrix is written to a file and can be evaluated by ERA or 
by a computer algebra system like Maple. 

The next section will give some examples of the identification situation in 


Factor Analysis Models. 
@ 
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4.4 Examples 


Consider the following model: 


YW = EArt Erot S3A3 +41, 
Yo = Sarqt E35 +€2, (4.4.1) 
¥3 = E3Ag +€3, 


where €,, 2 and €3 are unobserved variables. Instead we observe: 


xy = é, +€4, (4.4.2) 
XQ = &2 +és5, 
x3 = £3 +6, 


where the errors ¢,, i=1,..,6 are mutually independent and also independent of 
&1,-.,€3- All variables are assumed to be normally distributed. Let D be the 
covariance matrix of 4j,..,y3,%j,-.,%3, ® the covariance matrix of E15 ++5€3 
and W the diagonal covariance matrix of errors ¢,..,¢€s. The covariance 
structure is thus given by: 


Ay Az Ag vy 
As Orr Pi2 Pra 4,0 0100 v2 
yi i i A diz $22 $23 v2 Ag 0 010] 4+ Ys 
SSE 
di3 $23 P33 Az As Ag 001 va 
010 vs 
00 0 bs 
= ADA'+Y. (4.4.3) 


Here the number of equations, i.e. the number of subdiagonal elements in D 
(15), exceeds the number of free parameters (12) in A and 6: 6’ = (64 i. Pi3 
G22 $23 $33 Ar Az Az Aq As Ag). Consequently, the model might be identified. 

The Jacobian: matrix can be constructed using ERAFAC. The input-matrices 
that ERAFAC uses are the three parameter matrices above. Using these matrices 
as input, ERAFAC returns a 15x12 Jacobian matrix, where the first 8 columns 
are given by 
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0 Ag AAS ArAg AsAgtAdAs AsAS Pisrsthi2rAg Pa3sAsth22A4 











0 OAAs 29 Azvg AsX6 $136 $23A6 
Ar Az Az 0 ~0 0 P11 $12 
0 A O A Az 0° Pi2 $22. 
0 A, 0 Az 43 P13 $23 
0 0 0 0 joke tides 0 0 
© Ay Be 8 0 0 0 0 
0 0 de ‘4s 0 0 0 
0 0 0 0 he de 0 0 |, 
0 0 rA% 0 0 0 0 0 + 
0 0 0 0 Ay 0 0 0 
0 0 0 0 -_ 0 0 
0 1 0 0 0 0 0 0 
0 0 1 0 0 .0 0° 0 
0 0 0 0 1 0 0 0 
and the last 4 columns are "+ (4.4.4) 
bssrsthrar4 PosrstPorrothi2r1 Psarstbosrothi3A1 — 0 
“" b33r6 0 O $33Asth2srati3A1 
x3 - 0 0 0 
| Baa 0 0 0 
7 dss 0 0 0 
0 $23A6 b33r6 = PagAsthasA4 
* 07 Pri2 P13 0 
i 0 $22 $23 0 
i 0 $23 $33 0 
0 0 0 ; 13 
0 0 0+ $23 
0 0 0 $33 
0 0 0 
is, 0 0 0 
i 0 0 0 
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When evaluating this Jacobian by ERA, it gives a one-dimensional null-space of 
J(9), so that J(9) has a deficient column rank for all regular points, and 
consequently the model is underidentified by Corollary 4.4.1. The null—vector, 
n(@) say, being the tangent vector of the solution set {(¢., ¢,2 $13 $22 ¢23 
33 Ay Az Az Aq As Ag)}, can be presented as follows. 


- 1 go 2? 33 +¢ 1 $23 $23+92912033-291 2P13P23tPi3hisPo2 


oo oOo 0 © 


n(9) = | $22633A1-$23%23A1 . (4.4.5) 
—$12933A1tP1 39231 


$120 23A1-b1 39221 
0 


0 
0 


Even though the model as a whole is not identified, the parameters (¢1. ¢43 
$22 be3 $33 Aq As Ag) are identified. In order to ensure thé identification of 
the model as a whole, the first equation in (4.4.1) should have an additional 
restriction, since all parameters in the first equation are not identified. 
For example, the restriction A, = 0 would probably render the model 


identified. So consider, instead of model (4.4.1), the model 


MW = &Ar + §3A3 +€1, 
Yo = EaAgt E35 +€2, (4.4.6) 
¥3 = &3Ag +3. 


When computing the rank of the Jacobian matrix corresponding to this model 
(constructed by ERAFAC or by eliminating column 8 and equating A,=0 in 
(4.4.4)) one finds, indeed, that the model is identified. 
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; cd 
4.5 Conclusion 


Most of the earlier published results on restricted ‘factor analysis models 
pertain only to the uniqueness of A if W is assumed to be identified. For 
example, conditions for the local -uniqueness of A in orthogonal FA. (where © is 
restricted to equal J) were given by Dunn (1973) and Jennrich (1978), assuming 
the prior identification of ¥. Bekker (1989) gave conditions similar to the 
conditions in this chapter, treating the model as a whole. He also suggested 
how to proceed in evaluating the rank condition (see also chapter 3). 

Here a practical method has been presented that actually evaluates the 
identification situation of the parameters of restricted factor analysis 
models. The program ERAFAC takes care of the tedious job: of constructing a 
suitable Jacobian matrix and ERA computes its rank for regular points. 
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CHAPTER 5 


RESTRICTIONS ON THE REDUCED FORM: 
AUTOMATIC PARAMETERIZATION OF THE REDUCED FORM 


5.1 Introduction 


The precise number of constraints on the reduced form of a simultaneous 
equation system as implied by arbitrary exclusion restrictions on the 
structural form has been discussed by Bekker and Dijkstra (1990). In this 
chapter an algorithm due to Bekker (1986c) and a computer program will be 
presented that fully characterize these constraints on the reduced form. 


Consider a system of simultaneous equations: 
By = Tx +e, (5.1.1) 


where B is a non-singular mxm-matrix, I is an mxk-matrix and the disturbances 
€ are normally distributed and independent of x. The parameter matrices B and 
P are subject to exclusion restrictions, which can be presented in the 


following form: 
vec'{(B,I')}R = 0, " (5.1.2) 


where R is a matrix partitioned as 
R= , ; (5.1.3) 


and the matrices R;, ¢ = 1,...,m, consist of columns of the m+n-identity 
matrix in their natural ordering. It is assumed, without loss of generality, 
that R is of full column rank. Now let (B°,r°) be a matrix satisfying the 
exclusion restrictions (5.1.2), such that the ranks of all submatrices of 
(BP), also satisfying (5.1.2) are constant on an open neighborhood of 
(B°,r’). Then Bekker and Dijkstra (1990) prove the number of restrictions on 
the reduced form 7 = BP — whether or not the model is identified - to be 


equal to r,, where 
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7 = E in-p{(B°P RY], eS (5.1.4) 


p(.) denotes the rank, and n; is the number of columns of R,, t=1,....m. It 
should be noted that the number of restrictions on Mf is related to the 


manifold spanned by BT’ in the sense that 
Ty = mk — dim{(B"Y) | the restrictions given by (5.1.2)}. 


If the parameters of the i-th equation are identified, then p((B,P)R;) = m-1 
(see 3.2.3), and thus, if the parameters of all equations are identified, then 
r, equals Dj_,(n,-(m-1)) = Cyai n; — m(m-1). Consequently, in that case, the 
number of free parameters in JT equals mk—r, = m?4+ mk — Dij_, nj- m, which 
corresponds to the number of free parameters in (B, I’) after normalization. 

In section 5.2 a method for determining the nature of the restrictions on IT 
will be described. Section 5.3 will give an algorithm for the automatic 
computation of a minimal parameterization of IJ. This algorithm is used to 
create a computer program MINPARAM. Section 5.4 will give some examples using 
MINPARAM. Section 5.5 concludes. 


5.2 An algorithm for determining the nature of the restrictions on /] 


In this section the question of how the exclusion restrictions on (B,I’) 
translate into restrictions on J will be discussed. The answer will be 


formulated in terms of a minimal parameterization of II. 


Definition 5.2.1 
A local ‘parameterization of BP is a differentiable matrix function 118), 
with JT an mxk-matrix, 0 € R', so that BT and 41(@) describe, locally, the 


same manifold. 
Definition 5.2.2 


A local parameterization 1(6)=BT, Ge R', is said to be a minimal 


parameterization of BUD if Lis equal to the dimension of the manifold. 
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In order to formulate an algorithm for the characterization of the 
restrictions on J, we need another definition. Let y; denote the j-th column 


of P and let y;; denote the i-th element of y;. 


Definition 5.2.3 

Given an ordering of columns on (B,I), the restriction y;; = 0 is said to 
be transformable if $ is linearly dependent on the preceding columns in 
(Bo T° )R;. 
The number of transformable restrictions in the i-th row of If is equal to 
n,-p{(B°”°)R;}, so r, is equal to the total number of transformable 


restrictions. 


Example 5.2.1 


Let 
0 0 0 
By 0 0 0 V2 100 
B=} Bu B22 B23 |,F =| Ya 0 »R= ae , (5.2.1) 
0 fz Bag 0 0 000 
0 
then restriction y,,=0 is transformable, since | y2,| depends linearly on the 
0 0 x 
columns of | B22 623]. It can be checked that 73, = 0 is also a transformable 
B32 Bs 


restriction. These two restrictions are the only transformable restrictions 


that can be found in this model. 


Equation (5.1.4) shows that the restrictions on JT are the result of the linear 
dependencies that exist between the columns of (B’, i*, These linear 
dependencies are, in their turn, fully determined by the zero-restrictions on 
(B, I). The exact relation between the zero’s in (B, I) and the linear 
dependencies between the columns of (B’, r ) is given by the following Theorem 
that was proven by Bekker and Dijkstra (1990). 
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Theorem 5.2.1 Fe: 

Let 1; be the j-th column of r°. Let F° be a matrix consisting of different 
columns of ine Yooenss ¥j-1}- Then 1%; is linearly dependent on the columns of 
F° if and only if there exists a matrix F, consisting of n different columns 


of F° such that: 


(i) m-n rows of {Fh, 75} are zero, and 


(ii) no m-i rows of +4 columns of Fe are zero, i=0,...,n—1. 
If such an Fe exists, it has full column rank. 
Example 5.2.2 


In (5.2.1) we find that r is linearly dependent on the first column of B® 
and that of f°. That is, in terms of Theorem 5.2.1: 


0 
By, 0 12 
0 
i, = Bai Ya |: 12 =| 0 
0 0 0 


The third row of (Fe, 2) is equal to zero and no submatrix of Fe formed from 
m-—i rows and i+1 columns, i=0,1, is equal to zero. Thus the submatrix formed 
from the first two rows of Fe is non-singular and y is linearly dependent on 


the columns of a 


We can also use Theorem 5.2.1 to determine whether or not a particular 
restriction y,; = 0 is transformable. In section 5.2.3 an algorithm for a 
computer program will be presented that locates these transformable 
restrictions ‘automatically, and it also determines a matrix Fe, corresponding 
to this restriction as described in Theorem 5.2.1. 

Using the linear dependencies between the columns of (B’, a we will show 
below how a local minimal parameterization of JT can be found for points 
IT (0)=B TP in an open neighborhood of the regular point (B°, Ps. We start our 
procedure by parameterizing the first column of I: Bly,. Then we use this 


parameterization of By, to parameterize the next column of I: By. 





Together these columns provide a parameterization of eg” (V1, Y2). We proceed 
by using this parameterization to parameterize the next column of TT: B ys, 
and so on. In each step of this sequential procedure we may compute the number 


of additional parameters we need to parameterize the additional column by 
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using the expression given for r, in (5.1.4). We will show that the number of 
additional parameters to parameterize an additional column yy, say, is 
given by m, the total number of elements of BY; minus the number of 
transformable restrictions related to the corresponding column y; of P. 

Each transformable restriction Yiz=0, related to y,;, corresponds to a set 


of columns Fe, as defined in Theorem 5.2.1, so that 
Y= Fyp , 
where ? is a vector of size n. Now we may distinguish three cases: 


(i) There are no transformable restrictions related to Y;, in which case 
we may take Fe = B°. 


(ii) All transformable restrictions related to Y; can be characterized by 
exactly one set of columns | ol That is, for each transformable 


restriction y,;=0, the columns of Fe are also columns of (B®, PR, 
¥] a 


(iii) We need more than one set of columns Fe to characterize all 
transformable restrictions related to Y}: 


(i) 
If there are no transformable restrictions related to vy; then we have to use m 


additional parameters to parameterize Bohy, In other words 7; is not 


restricted: every element is a parameter. 
(ii) 
In this case we have 
15 = Fr’, (5.2.2) 


where n, the number of elements of p°, is simply m minus the number of 


transformable restrictions related to Yy 
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Now assume that Sg and Sp are selections of. columns of Im and J, such that 
Fe = (B°Ss, rsp). Since (B’, r°) is a regular point, it- can be stated that, 
for the points (B, I) in an open neighborhood of (B°, ar ), 


p(BSp, TSr) = p(B°Sz, I'°Sr), (5.2.3) 
and 
p(BSp, Sr, 1;)= p(B°Sp, T°Sr, 75)- (5.2.4) 


This implies that the linear dependency (5.2.2) also holds in an open 
neighborhood of (B’, ry, That is, for every point in this neighborhood, there 
exists a vector p so that: : 

Yj = (BSpz, TS y)p. (5.2.5) 
Premultiplication of (5.2.5) with B' gives 


n; = (Sp, MSp)p, (5.2.6) 


which shows that the j-th column of 7 can be parameterized using the 
parameters in the preceding columns of J and n additional parameters in p. 


Example 5.2.3 ; 
We may now construct a minimal parameterization of the matrix I = Br’, 


where (B, I) is restricted as in (5.2.1), . ase @ 
By 9 0 0 12 j ae 
B=} Ba B22 Basl|, =| Yn 9 |- an (5.2.1) 
0 B32 Bs 0 60 


This model has two transformable restrictions: y,,;=0 and ‘Y23=0. The 
transformable restriction y,,;=0 corresponds to the linear ‘dependency: 
€ . . 


0 60 


%1 = | Boz Bos | Pr acs ; (5.2.7) 
B32 B33 
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Let p; = (6,,6,), then 
0 0 0 
6 
™m={1 0 (e] = | 6, |. ; (5.2.8) 
01 8, 


The transformable restriction 3.=0 corresponds to 


By 0 


Yo = | Ba Yar | Po- (5.2.9) 
0 0 


Let py = (63,04), then 


1 0 8s 
t= | 0 8, (@°| = | 0,8, |. (5.2.10) 
0 6, 0.04 


Thus a minimal parameterization of II in an open neighborhood of (B®, py will 


be given by 
(8) = | 6, 6,04}. (5.2.11) 
8, 0264 


Indeed we find that the number of parameters used in the minimal 


parameterization equals mk, the number of elements in I, minus 7, = 2. 
(iii) 


Now consider the minimal parameterization of a column m; where more than one 
set of transformable restrictions is related to 73, that is, where we need 
more than one set of columns ye to characterize all transformable 
restrictions. If the linear dependencies between 1 and the preceding columns 
in (B°, ry can only be characterized by more than one matrix F,, say Frys 
Frio where Fa,= (BSz,, Sr ,); i = 1,...,r, so that 7} = Fa Teim % = a at 


then we may state: 
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Lr (Sz, MSy,)P,> 


(5.2.12) 
m; = (Sp, Sp, )p_. 


So, 7; can be found as the intersection of (Sz,, MSr,); t=1,...,7. The basis 
of this intersection can be parameterized with parameters that have been used 
to parameterize the columns 7,, t<j, while the dimension of the intersection 
equals m minus the total number of transformable restrictions related to 7;. 

A general method to compute a parameterization of the relevant intersection 
is as follows. Let N(A) be a parameterization (with the parameters of A) of a 
basis of the null-space of A, so that AN(A) = 0. Then a local parameterization © 
of a basis of the intersection of A, and A, can be found by (A,,0) N(Aj,A2). 
In fact, the parameterization of N(A) is exactly what ERA (see chapter 3) 


computes: the null-space of a parameter—matrix. 


Example 5.2.4 
Let the restrictions on B and I be given by 


ete (PLE TL ge 

B= 21 P22 Pas pes Y22 ‘ (5.2.13 

Bs2 B33 Bsa!’ 0 0 Y33 0 ( 
0 O Baz Bas 0 0 ODO ¥s4 


With respect to the third column of [: 7,3=0 and 4;=0 are transformable 


restrictions, so: 


0 0 0 
_ | B23 0 Y22 of 

Te B33 Bsq 0 iy : : Pei) 

Bas Bag 0 
and : * 

Bir Bia so 

ii fa Sa ; Bs . (5.2.15) 
0 0 0 


As a result of (i) and (ii) the first and second column of HT can be 


parameterized as: 
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9, 9206 
8, 9396|" 
84 4496 


(5.2.16) 


Then, 


1; = (5.2.17) 


oroo 

rPooce 
iv} 
ko =) 
o 


and 


T; = Pay (5.2.18) 


oocr 
coore 
S 
2) 


It follows that 


Os 0 
0,6, 0 
tt, = a C 05 p. (5.2.19) 


Indeed we find that we need m minus the two transformable restrictions in the 
third column of I, i.e. 4-2=2, additional parameters to parameterize the third 


column of JI. Let p’=(67,0g), then the third column can be parameterized by 


050, 
x, = | 92969 
2= | 6565 


9498 


(5.2.20) 
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5.3 An algorithm for the computation of a minimal parameterization of /I 


In this section it will be shown how the theory in section 5.2 can be 
implemented in a computer program. The next section will give some examples 
using the computer program based on the algorithms given in this section. 
Schematically a computer program that returns a minimal parameterization of I, 


given the restrictions on B and I’, may look like this: 


1. Read the restrictions on B and [ 

2. Test for all restrictions y;;=0 whether they are transformable. 
If a restriction is transformable, save a corresponding set eS. 

3. Construct JT using the information saved in step 2. , 

4. Write I. . 


Of course each step needs to be refined. Here we will only refine step 2, , 
since it is the most complicated step. How step 1, 3 and 4 should be refined 
will not be discussed here: steps 1 and 4 are trivial, and step 3 has already 
been discussed in the previous section. Step 2 can be refined as indicated in 
algorithm 5.3.1. The table provides some definitions of variables used in this 
algorithm, which might prove helpful in understanding the algorithm. 


Table 5.3.1 


variable function/explanation 


j current column number of I’. 

a current row of (B, I), i.e. the th equation. ; 

F’ is a different set for each i, it is a selection of columns 
from (B, I)R;. The initial value of F * is a selection of 
different columns of B, depending on R;. In the j-th ‘step it is 
a selection of columns of (B, 7j,...,7;-1), depending on R,. 





m total number of rows of B and I’. 
é k total number of columns of I. 
n current set size used to locate a set F,,: . 
Fy, current subset of n columns of F* 
, Sy will hold all solutions (F,,,7;) in terms-of Theorem 5.2.1 
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Algorithm 5.3.1 


21 


2.2 
2.3 
2.4 
2.5 
2.6 


2.7 


2.8 


2.9 


2.10 
2.11 


[Initialize] 7 <— 0, rz <— 0, solution set S; < {}. 

Create m sets F i which contain the columns of B, that also appear 
in (B,P)R;. 

[Loop on j] i <0. j —j+ 1. if j > k, then go to step 2.9. 

{Loop on i] i <— i+ 1. If i > m, then go to step 2.2. 

[Exclusion restriction?] If y,; # 0, then go to step 2.3. 

[Initialize] n < 0. 

[Loop on n] n <n + 1. If n > number of columns in F. then y,; =0 is 
not transformable: go to step 2.10. 

{Construct F,,] Let F, be a subset of n columns of F *. Create an F,, 
that has not been tested by step 2.8. If such an F,, does not exist 
(ie. all possible subsets of size n have been tested), then go to 
step 2.6. 

[Test dependency] If y,; is not linearly dependent on F,, then go to 
step 2.7. 

[Solution] ,; is transformable. We have found a subset Fe = F, for y; 
as in Theorem 5.2.1. Add the solution (F,, Yj) in the solution set Sy. 
Tx <— Tx + 1. Go to step 2.10. 

[Increase F*] F < F* + {y;}. 

[Finish] 7, contains the number of restrictions on IJ, Sy contains all 


solutions. 


Steps 2.7 and 2.8 should be refined even further. For example, the program 


could use a recursive procedure to construct all possible subsets of size n 


from F . and step 2.8 could be programmed using Theorem 5.2.1. However, these 


two steps are very inefficient: for each constructed set F,, all the necessary 


conditions of Theorem 5.2.1 have to be checked. A more efficient method is 
based on Theorem 5.3.1. Let X be a r(X) x c(X) matrix of rank p(X), with z(X) 


zero rows and z(X) = r(X) — z(X) non-zero rows. 


Theorem 5.3.1 
Let F° be a matrix of full-column rank, consisting of different columns of 


{B°, Pagveey 75-1}. Then Y; is linearly dependent on the columns of F° if and 


only if there exists a submatrix F. of F° so that 2 (Fe, 75) =n. 
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Proof 


(=) If 2 (Fasvj)=n and Fe is of full column rank, then without loss of 


generality, we may write in partitioned form: 
(Fe, 75) = | Pg } with P a non-singular nxn—matrix. - 
0 0 


It follows that 75 is linearly dependent on the columns of F,. 


(=) If 3 is linearly dependent on the columns of Fe then according to 
Theorem 5.2.1, Fe exists so that z(Frs¥3) = m-n, and thus z (Fe, a =n. 


We can improve the efficiency of the algorithm dramatically by employing the 
consequences of Theorem 5.3.1. 

First, ensure that F’ is such that any matrix formed from n columns of F ; 
has rank n. This can be done by replacing step 2.9 by: _ 
2.9’ [Solution] ,;; is transformable. Add the solution (Fn, Y;) in the 

solution set Sy. Tz <— Tz + 1. Go to step 2.3. 


As a result ; will not be added to the set F oe Yij = 0 is a transformable 
restriction. Indeed this enlargement is not necessary. Let X be a subset of 
columns {7j41,---;Yt-1} and let y; be linearly dependent on the columns of F* 
(Yij is transformable). If a restriction iz, t>j, is transformable, then 7; 
is linearly dependent on the columns of F 7 {y;} U X. So 7; is also linearly 
dependent on the columns of F Deh X, which is a smaller set. As a consequence 
we know that any matrix F,, formed from n columns of F* has rank n, n<xm. This 
is also true for the initial values of F* , Since B is a non-singular matrix. 
Secondly, we know that the minimum number of columns n in F,, is equal to 


2 (7;)- Consequently, step 2.5 can be changed in: 
2.5’ [Initialize] n <— z(y;)-1. 


Third, since we know that any matrix F,, formed from n columns of F* has 
rank n, n<m, step 2.7 and 2.8 can be changed in: "find a matrix F,, such that 
Z(Fp, Yj) = n”. Furthermore, if z(F,, y;) =” and F; consists of a selection 
of t different columns of F,, then z(F;, 7;) <n. Hence we can eliminate in an 


early stage all sets F,, containing F,, with n < Z(F,, 13) 
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Let F’ be a selection of columns of (By Yi, ++) Yj) as in algorith: 
5.3.1. Let this set consist of, say, c columns {f,,f2,..,f-}. Now a subset F, 
of size n has to be found such that 7; is linearly dependent on the columns of 
F,. The algorithm below describes how such a F,, can be found. First this 
algorithm will be described in words. 

For a fixed n, F, can be formed in a recursive stepwise procedure. Let, in 
a recursive procedure, F * be a subset of F'. Initially we take F * to be empty. 
Then in each step we consider an additional column, f; say, of F * not 
contained in F’. If Z(F , fis Y;) <n, then f; is added to F’, otherwise 
another column will be considered. If no f; can be added, then the last column 
added to F’ will be deleted from F’, and the procedure is repeated. In fact, 
no more columns f; will be considered as soon as the number of columns in F’ 
that remain to be considered, is less than n-c(F’). If finally Z(F’, fi 13) = 
n, then F F F,, and y;; = 0 is transformable. Otherwise —- we have considered 


all (*] combinations - we fix n to another value: n<—n+1 (step 2.6). 


Now a more formal description of the algorithm follows. We use the variable 
names as in table 5.3.1. Additionally, variable t¢ indicates the number of 
columns currently in Fs and F” a subset of F’, as indicated above. 


Algorithm 5.3.2 


2.7.1 [Initialize] F, — {}, F’ <— {}, t — 0, 1 <— 0, stack < {}. 

2.7.2 [Loop on t] | <1 + 1. If lin-t > c+1, then go to step 2.7.7. 

2.7.3 [Test] If 2(F, fi. Y;) > ”, then go to step 2.7.2. 

2.7.4 [Add column] F’ < F'+ {f,}, t — t41. 

2.7.5 [Recurs] If ¢ < n, then push / on stack and go to 2.7.2. 

2.7.6 [Solution: t=n] F, <— F ie clear stack and finish. 

2.7.7 [Recurs back] If stack is empty, then finish, otherwise pop | from 
stack, FP an Fo {fi}, ¢ <— t-1 and go to step 2.7.2. 

2.8 [No solution] If F,, = {}, then go to step 2.6. 


The following example provides an illustration of this algorithm. 


Example 5.3.1 


Fe= {fi, fa fs fa}, n = 3. 
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xy X4 X7 Xs X10 
0 X5 0 X9 
fi= | ol, fo= | el, fs=1]9),f4=] 01,25 = | Xu 
0 0 0 
X3 0 0 0 0 


In each line the step of the algorithm above and the corresponding values of 
the different variables will be given. A dot indicates that a value is 


unchanged in the current step, so the previous value still holds. 


Step t t n F Z(F ,fyy;) stack - Fp 
2.7.1 0 0 38 {} 3 {} " {} 
27.2 1... . ; ; ; 

2.7.3 ee ks . 3 

2.7.4 1. {fi} : 

2.7.5 . 4. . ; {1} 

2.7.2 2 : : 

2.7.3. 4 : 

2.7.2 3 ; ; t. 

2.7.3 — 4 2 ’ 

2.7.4 2. {fifa} ‘ : . 

ok nn aor , ’ {1,3}... ‘ 
2.7.2 4 : . é 
2.7.3. : 4 3 

27.2 5 . . , ; ; 

27.8 3 1~. {fi} : Fale! 

27.2 4 .. 4 : . 

2.7.8 1 0. {} : {} 

27.2 2 ... ; i . 

27.3 . 2. . 3 

2.7.4 iy: {fo} 2 

27.5 . . . : {2} 

2.7.2 3 ; 

2.7.3 es ; 3 ; 

2.7.4 , 2 =. {fo fs} : 

2.7.5 ee : {2,3} 

27.2 4 ... : : ‘ s 
273 = es ' 3 ; e 

2.7.4 3. {fo fa,f a} : : 

2.7.6 A . : {} {fafsrfa} 


So y; is dependent on {f2,f3,f4}- 


Now our algorithm is complete. A computer program - MINPARAM - has been 


created using the algorithm as indicated above. 
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5.4 Examples 


Using the computer program MINPARAM we are able to construct some examples of 


a minimal parameterization of 77. 
Example 5.4.1 


Consider the model where 


Bu Bio - : Yn - 
Ba Bo2 Box - » Yes « 
- B32 Bs3 Bq ; #5 coe) Beare IP (5.2.13) 
Bas Baa im 2 Fag 


Given the restrictions on these matrices as input, MINPARAM returns as output 
the following information: 


(i) The number of restrictions on If: rz, = 6; 


(ii) The linear dependencies between columns of r° and B’, indicated as 
below, where each row specifies a dependent set of (B°,r’), 


Yo Br oY 
Ys Bs Ba Yu 
Ys Bi B2 Yu 
Ys Bs Bs Yn 
Ya Ba 3; 


(iii) A minimal parameterization of the matrix I: 


9, 95 9567 956789 
0, 9296 929687 8296978 
83 8386 930, 939885 
94 9486 9483 P10 


(5.4.1) 


If the diagonal elements of B are normalized then the number of unrestricted 


parameters in (B, [) equals the number of parameters in //: 10. This indicates 
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that the model is identified after a normalization. 
In this relatively simple example we can show in another way that this 
parameterization is correct: by using a computer algebra system, like Maple, 


we can compute B'r. The resulting output is given by 
I = A” P, ; (5.4.2) 
with 


A= 8118 228338 44-F 118 229348 43-11 328 239 44-P 21 128 338 44+ P2112 348 43 


and where the first two columns of P are parameterized by, 


(B28 338 44-8 22P 34843-32823 44) V1 —Pr2( 8338.44—P34P 43) 22 


—Ba1 (833844—-834P43)111 811 ( 8338 44—Bs4P43)¥22 
BoB 328 44711 —B118 328 44722 : : 
—Bo18 328 43V11 B18 328 43722 


the last two columns of P are given by 


8128 238 44733 ~Bi2h 238 34744 
-B118 238 44733 B18 238 347 44 

Ba4( 8118 22-BarPi2)¥33 —Bsa( B11 822-B21Pi2)¥44 { 
—B43(B11822—-B2rPi2)¥33 (8118 228 33-8 118 23832—B 218 128 33) ¥44 


Equating this parameterization with (5.4.1) shows that 


9, = (228338 44-P22834P3-Ps28o3h aa) / A, 
92 = —B2;(B33844-PsaBaz)¥u / A, 

93 = Bar Bs2Batu / A, 

94 = —Bahs2Bax¥u / A; 

95 = —By2(B338 44-P34Pa3)¥22 / 4, 

95 = -Bul Bir, 

8; = —B32B44/ (8338 44-P 3484s), 

63 = — (By 822-B2i812) / Bubs, 

9 = -Ba3/ Baas : 
19 = (8118 22833-B 118 238 32-219 12833) 44 [A . 
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This shows that the minimal parameterization is indeed correct. On the other 
hand we find that the minimal parameterization (5.4.1) is much less 
complicated than the parameterization based on (5.4.2). Moreover, the latter 
parameterization will be minimal only if the model is identified. 


Example 5.4.2 


a. Let the model be specified like: 


Bu. Bis... Bay Yn. 
B22 B23 Bos Bas Bre . foe (> ee 
B33 B34 83s B36 Bsr ee WEB mes 
B= |. . Bas Bas Bas Bas Barii T= 1. 2 | Meet; (5.4.3) 
Bss Bse Bsz c « » & Fe, ; 
Bes Bes - sw os ow a) OS N67 
Brs » Baz Pa ree 


The output of MINPARAM is given by: 


(i) The number of restrictions on HI: r, = 22; 


(ii) The linear dependencies between columns of °° and B°: 
1 Pry 
Yo Ba 
Ys Pi Bo Bs Bay 
Ys Bo Ba Ys, 
Ye Bo Bs Bs Be Y¥3 Ys 
Y7 Be Bs Bs Be Ys Ys, 
Yr Bi Br Ys Ya Ys Yo 
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(iii) A minimal parameterization of the matrix Z : 
columns 1 to 6: 


9, 0 83 8367 O15 93017+810918 
8, 94 9, A Aro 


0 
0 0 85 A587 812 956174912918 
0 0 06 9 3 Bao 
0 0 0°90 614 924 
0 60 615 9a 
0 0 M16 4618 
column 7: (5.4.4) 


93679 2 4+9 39274939 59) 7823-83828 8923+ 93012925+93912926+959106 25 
0482746 59; 192640591 9923+ O24 

95828 » 4+ 958974 58 5917923 +059129 1892395912926 

9581382 6+85929923+96927+8 824 

959149 26+95821 923 

9581582 6+85922903 

9:616925 


We find that IT is parameterized with 27 parameters, which equals the number of 
free parameters in B and I (after normalization). Thus, the model is, again, 
identified. ; 


b. On the other hand, let the model be specified like: 


Bu. Bis... Bar TAs 5 
B22 Bos Baa Bos Poe . Hh ke ee 
B33 B34 B3s B36 Baz A ee) a | 
B= |. . 8a3 Bas Bas Bas Baris T= |) 2 2. Yue Ok 
Bss Bse Bsz Se ee ME 
Bes Bes Bex > hk ae . Yes. 
Bzs Bre Brz me a ae Oe 
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We find the following MINPARAM-output: 
(i) The number of restrictions on I: r, = 22; 


(ii) The linear dependencies between columns of °° and B°: 
n By 
Yo Ba, 
Ys Bi Bo Bs Ba, 
Ys Bo Ba Ys; 
Ye Bo Bs Bs Be Ys Ys; 
Yr Bo Bs Bs Be Ys Ys 
Yr By Br Ys Ya Ys Yoi 


(iii) The same minimal parameterization as in (5.4.4). 


So this model implies the same restrictions on IT as the model in a. We thus 
find that the models (5.4.3) and (5.4.5) imply the same restrictions on JZ. In 
that case the models are said to be (locally) equivalent (cf. Bekker, 1990): 
the data cannot be used to distinguish between the two models. 


c. Consider model (5.4.3) where the last column of P has been deleted: 


Bu. Bis . ‘ . Br Yu. : 
B22 Bos Boa Bas Bre . e 222 5 ‘ 
833 B34 B35 B36 Bsr s » 38 5 : ‘ 
B= |. . Bas Baa Bas Bag Bar|3 T= | . : . Yaa. : (5.4.6) 
8s5 Bse Bsr ie 4s ‘ ; VSS « 
Bes Bes - % ; ‘ . 66 
Brs - Braz % = ‘ ‘ . 76 


A minimal parameterization of JT is simply found by deleting the last column of 
Mf in (5.4.4). So Hf is parameterized using only 22 parameters. Since in this 
model there are 25 parameters to be estimated (after normalization), the model 
is not identified. However, we know all restrictions on JT and as a consequence 
IT can be estimated efficiently, even though the parameters of B and I are not 
identifiable. 
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5.5 Conclusion 


In this chapter it was shown how the computer program MINPARAM constructs a 
minimal parameterization of IT. This parameterization is found by scanning the 
linear dependencies between the columns (B, I), in case (B, [) is subject to 
zero-restrictions. This minimal parameterization can be used for the efficient 
estimation of the reduced form parameters whether or not the structural form 
parameters are identified. 

The next chapter considers general linear structural models, where we will 


find further use for the minimal parameterization of matrices like ar 


90 








CHAPTER 6 


COMPUTERIZED EVALUATION OF RANK CONDITIONS: 
IDENTIFICATION IN LISREL MODELS 


6.1 Introduction 


In chapters 3 and 4 we presented methods to evaluate the identification 
situation of simultaneous equation models with measurement error and factor 
analysis models, respectively. These models are both subsets of the so called 
LISREL model (Jéreskog, 1973). In this chapter we present theoretical and 
practical methods to characterize the identification situation in LISREL 
models. 

Seidel (1988) discusses a computer program — LISRAN —- which constructs and 
evaluates a Jacobian matrix to characterize the identification situation of 
general LISREL models. The method he used for the construction of the Jacobian 
matrix is rather straightforward: the derivative of the symbolic covariance 
matrix is taken with respect to the free parameters in the model. LISRAN then 
evaluates the rank of this symbolic matrix in order to characterize the 
identification situation of the LISREL model. However, the Jacobian matrices 
constructed by LISRAN tend to consist of elements containing many additive and 
multiplicative terms. This makes these symbolic matrices not very suitable for 
manipulation by computer. 

In this chapter we will derive a number of Jacobian matrices, each suited 
for evaluation of the identification situation in LISREL models. Some of these 
Jacobian matrices can only be used if the LISREL model under consideration has 
specific a priori restrictions. We also present computer programs that 
construct these symbolic Jacobian matrices. The ranks of these matrices can be 
evaluated using ERA (cf. chapter 3), or by a CAS like Maple. 

In section 6.2 we present the general LISREL model and show that it can be 
transformed into a restricted FA model. In section 6.3 we use a minimal 
parameterization in order to characterize the identification situation of a 
particular LISREL model. In section 6.4 and section 6.5 we derive Jacobian 
matrices that can be used to characterize the identification situation in 
general LISREL models. The Jacobian matrix given in section 6.5 does not have 
any additive or multiplicative elements; it has, in our view, the most 
suitable form for evaluation by a computer algebra system. Section 6.6 


concludes. 
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6.2 The LISREL model considered as an FA model . | 
il 


Consider random vectors 7’ = (5Ma5-5%m) and €'= (€162)-y€n) of ldtent 


dependent and independent variables, respectively, and the system of linear 


structural equations 


Bn + TE = 6, (6.2.1) 


where B is a non-singular mxm-matrix, [ is an mxn-matrix and ( = 


(C1,Cay+9Gm) @ random vector of residuals. The vectors » and € are 


unobserved; instead we observe y' = (¥1,Va-s¥p) and x'= (%1,%,..,%q), such 
that: 
y=A,nt+é, (6.2.2) eae 
x=h E+ 6, (6.2.3) 


where A, is a pxm-matrix and A, is an gxn—matrix. It is assumed that 
1. , €, € and 6 are mutually uncorrelated, 
2. B contains “ones” on the diagonal and is non-singular, 


3. All variables are measured in deviation of their means. 


Let Dg (nxn) and - (mxm) denote the covariance matrices of € and ¢ 


respectively and ©, and $s the covariance matrices of € and 6. 


Formula (6.2.1) can be rewritten as 


n= (B? ,-B'D) | i | (6.2.4) 
is 





and (6.2.2) and (6.2.3) can be reformulated as 


P| ag | ay OTP BP |e) a) ed (6.2.5) 
x 0 AJLO In é 5 
or, 
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Ber Ge 


Now let 


an(4oha=[? 15 [2] t= [£] me [502m 
0 A, 0 Ip x é 5 


then (6.2.6) can be rewritten as 
y= AA Ese. (6.2.8) 


If the parameter matrices B, IT, A, and A, are restricted only by zero 
restrictions (and normalizations), then we can use the method described in 
chapter 5. That is, the computer program MINPARAM can be used in that case to 
compute a minimal parameterization of the matrix 


A= AA”. (6.2.9) 
Using such a minimal parameterization we can further reduce (6.2.8) to 
yu Aése, (6.2.10) 


which is an FA model with general restrictions on the parameters; in 
particular, A contains non-linear restrictions. 

The identification situation of such general FA models can be evaluated 
using the methods of chapter 4. We can use the method (ERAFAC) described in 
that chapter to investigate the identification situation of the transformed 
LISREL model. If the transformed model (6.2.10) is identified according to 
ERAFAC, and the total number of parameters in the minimal parameterization of 
A equals the total number of parameters in A and A, then, as can be easily 
proved, the original LISREL model is identified as well. 

So, the identification problem of a subclass of LISREL models — where the 
matrices B, I’, A,, A, are subject to zero restrictions -— can be handled, 
without any modification, by methods already described in chapters 4 and 5. In 


the next section we will show that a minimal parameterization could also be 
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used to characterize the identification situation of LISREL models which 
contain only zero restrictions on B and A,, and general restrictions on Py, A, 


and the covariance matrices. . 


6.3 Identification of LISREL models using a minimal parameterization 


In order to describe another approach for characterizing the identification 


situation of LISREL models, we rewrite formula (6.2.5) as 


seeded 


The second-order moment equations are given by 
XZ = AANA'A' + &, (6.3.2) 


where ¥ (ix), l=p+q, Q (kxk), k=m+n, and = (ixl) denote the covariance 
matrices of (y',x’), (¢’,é’) and (e¢’,6’), resp., and 








A=l|P 9} anlie P| pel % 9} anaza| % © | 633) 
0 A e x 0 % 0 8; 
with P = A,B". 


If the restrictions on the elements of A, and B are of the exclusion type 
only, then we can find a minimal parameterization of Bd, by using the 
computer program -— MINPARAM — given in chapter 5. This minimal 
parameterization may give a first indication of the identification of the 
LISREL model: it shows whether or not the parameters in B and Ay are uniquely 
determined by the parameters in P. If these parameters are not uniquely 
determined by P, then they cannot be identified in ‘the LISREL model. However, 
if the parameters in B and A, are identified by P, then the identification 
problem in LISREL is reduced to the identification of the parameters in A, A, 
2 and P. 

Now we can take the derivatives of Z with respect to A, A, 2 and =, 





8 vec Y 
8 vec'A 





= (T2+K;,,1) (AARA'eI)) > (6.3.4) 
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8 vec 


= (F2+K;,1) (AeA) (24el;) > (6.3.5. 
8 vec’A 
vec ¥ _ (hed) (A@A), (6.3.6) 
8 vec’N 
ova a Fe, (6.3.7) 
8 vec’= 


where K,; is a commutation matrix. The expressions (6.3.4)-(6.3.7) each 
contain (? rows, and kl, kK’, K’, P columns, respectively. 

The parameter matrices A, $, and = are linearly restricted, but A may 
contain non-linear restrictions due to the minimal parameterization of P. Let 
P=P(8), with @ the vector of minimal” parameters, and A,=A,(A), with A the 


vector of free parameters in A,. In that case the matrix 


8 vec A 


A ee NEC UE 
Ovec’(8’,r') 


is easily computed. Let the n,, say, linear restrictions on A, 2 and = be 
given by 


R'vec(A, 2, 2) = 1, 


where R’ is a constant matrix of order n,x(2k? 417), and r is a constant vector 


of order n,xl. The Jacobian matrix is now given by 


avec LY 
O(8’ A’ vec’(A, 2, F)) (6.3.8) 
a 8 R'vec (A, 2, =) = 
0,0, ——_ 
Ovec'(A, 2, =) 
(1y2+K;,) (AARABI,) A, —(T2+Ki,1) (AA2@A) , (A@A) (A@A), T2 
(6.3.9) 


0 : R’ 


Formula (6.3.9) could be programmed at face value. However, programming 
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Kronecker products and commutation matrices efficiently is rather tedious. It 
is much easier to create a computer program, for taking the derivatives of Y 
directly. So, let 9* be the vector of all free parameters in (6.3:3), then the 
computer takes the derivative of Y with respect to 9": 


_ Ovec ¥ (6.3.10) 


Jy 
r,) Q* 


Algorithm 6.3.1 provides a method for constructing and evaluating J. 
Algorithm 6.3.1 Construction and evaluation of J. 


Read parameter matrices A,, A,, 2, B, f and &. 

Construct output of B’, A,, suitable as input for MINPARAM. 

Call MINPARAM with the input created in step 2. 

MINPARAM returns the parameterization of P’= B’*A,. ; 
The output of MINPARAM can be used to see if the parameters of B and Ay 
are uniquely determined by P. If the parameters of B and A, are not 
identified by P, then stop: the model is not identified. 


eek NS 


5. Read the parameter matrix P. 

6. Construct Y as in (6.3.2). 

7. Determine the vector 9, say, of all free (=different) parameters in the 
parameter matrices A,, 2, F and = and P. 

8. Take the derivative of © with respect to 6, and write the resulting 
Jacobian matrix to file. Write the contents of @ to another file. 


9. Call ERA (or Maple) with the Jacobian—output in step 8 as input. 
A computer program — LISPARAM — has been written based on algorithm 6.3.1 to 
evaluate the identification situation. An example of the characterization of 
the identification situation of a LISREL model using this computer program is 
given below. : 


Example 6.3.1 


Let a model be specified as 


1 Bi2 ny Yii Yi2 éi C1 
4: - ; (6.3.11) 
Ba 1 Ne 0 ol 2 C2 
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VW A, 0 Ey 
ya Az 0 ™m &2 
yz | =| Ag 0 + | z, |. (6.3.12) 
Ya 0 4 Ne &4 
¥5 0 As Es 
xy 1 0 6; 
& 
x. | =| Ae 0 + | & |, (6.3.13) 
Xs 0 1 &2 53 


with the variance of 6, = 0, ie. 6; = 0 almost surely. The matrices Gr, , 
and @5 are assumed to be diagonal, while $, is only symmetrical. Thus we have 
22 unrestricted but unknown (free) parameters. 

Following algorithm 6.3.1, we first build a minimal parameterization of 
B'A,. The computer program MINPARAM gives 


9, 8,6, 0:04 95 050, 
0, 9203 204 96 667 (6.3.14) 


as a minimal parameterization of P’. The minimal parameterization contains 
seven parameters, which equals the number of parameters in B and Ay. So, the 
parameters in B and A, are uniquely determined by P. Consequently, B and Ay 
will be identified if and only if the parameters of P are identified. 

Subsequently, LISPARAM takes the derivatives from £ with respect to {Yu, 
12, Ae, 91, G2, 83, 4, O5, O6, 47, @e os a £2.43 Pe os Peas o..,, 2... 
o,., o,,, $.., 5, 5}; which yields a 36 x 22 Jacobian matrix. 

When evaluating this Jacobian matrix by ERA we get a 22 x 2 null-space, 
indicated by 


oo 
oo 
oo 

*¥O 
oo 
oo 
oo 


(6.3.15) 


oo 
oo 
oo 


where a ’*’ indicates a non-zero value. Consequently this LISREL model is not 
identified. However the output also shows that not all individual parameters 
& & & 


as Und Dk DS 


are underidentified: {Ag 43, 94 4, Pe. %e,.., Pe & 
$,. Ds 


£59 %5,} are identified. Unfortunately we do not have such detailed 


ak 
information with regard to the individual parameters in B and Ay. Since P is 
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not completely identified we cannot indicate which parameters in B and A, are 
identified and which are not. : 

In his doctoral thesis Plat (1988) describes a model like (6.3.11)- 
(6.3.13), with the additional restrictions 


« 


A= Ay = 1, (6.3.16) 


Since the restrictions in (6.3.16) are non-zero restrictions, MINPARAM cannot 
be used to find a minimal parameterization of B’A, in this case. However, by 


hand we can show that the matrix 


1 0, 03 4 0495 
0, 0:0, 0,03 1 4s 


provides a minimal parameterization of B’A,. 

Using this parameterization we can start algorithm 6.3.1 at step 5. If the 
Jacobian matrix, constructed by the rest of the algorithm, is evaluated by 
ERA, no null-space is returned, and thus the model (6.3.11)-(6.3.14) with the 
additional restrictions (6.3.16) is identified. This was also concluded by 
Plat in the appendix of his thesis (1988, 130-135), using manual elaborations. 


6.4 Identification of LISREL models in the general case 


In the previous section we presented a method to investigate the 
identification situation in a subclass of LISREL models, where the 
restrictions on the elements of A, and B are of the exclusion type only. In 
this section, and in the following section, we will give methods to 
investigate the identification situation of general LISREL models, in the 
sense that rather general linear restrictions are imposed on the parameters of 
the model. 
Using (6.2.6) we can write the second-order moment equations 


D=AA'NA AS SE, (6.4.1) 





where ©, 2 and = are the covariance matrices of (y’,x’), (¢’,€') and (€',6’), 


resp., and 
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fe | 2) ae |PP) gl Se & (6.4.2) 
0 A 6 ¥. 0 &, 


Note that the definition of A and A in this section differs from the 
definition (6.3.3) in the previous section. 
Let the n,, say, linear restrictions on the parameter matrices A, A, 2 and 


= be given by 





i R’ vec(A, A, 2, =) = 1, (6.4.3) 


where R’ is a constant matrix of order n,x(ki +2k? +17), and r is a constant 


vector of order n,xl. The Jacobian matrix, 


d vec ¥ 

Ovec'(A, A, 2, F) 
@ R'vec (A, A, 2, F) 
Ovec'(A, A, 2, =) 


(6.4.4) 


should be of full column rank in order for the parameters to be identified. 
Using standard matrix calculus, we find the following expressions: 


VEC 
(2 + Ki) (A @ 1) (A72Ae 0), (6.4.5) 
avec Z % K,,) (AA72 @ AV (A™9 AR 
Ovec’A at ad wd C = ) ¢ ® ) = 
-(Ip + Ky) (A @ A) (4724 A, (6.4.6) 
ig = (A@ A) (479 AY, maa 
VEC 
dvec Y _ ie. es 
Ovec’= 


So, the Jacobian matrix is given by 
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(I+K} (Ae) (A‘2A™ eI), -(I+K) (AoA) (A*2A"*@A"), (AA) (A7*@A) I 
R’ 


(6.4.9) 

with K = Kj; the proper dimensions of J can be found in (6.4.5)-(6.4.8). 
This matrix is not very well-suited for evaluation since inverse matrices 
tend to consist of many additive and multiplicative terms. However, J can be 
simplified somewhat by postmultiplying it by the non-singular matrix A, where 


(A'n" Agel) 0 0 0 
0 1971 ; 
A (A’'R AeA) . 0 0 (6.4.10) 
0 0 AeA 0 
0 0 0 I 12 


which yields 


(11 2+K, , 1) (Aol) , ~(N12+Ki,1) (A@A) , (A@A), Tz 


A'2”*Agl, 0 0 
JA= - ; AG" AeA oo | | (64-11) 
0 0 AeA 0 
0 0 0 Tz 


The matrices J and JA both have Pin, rows and kl+2k’+l? columns. Since A is a 
non-singular matrix, the rank of JA equals the rank of J, and thus the 
identification situation of a LISREL model can be investigated by evaluating 
the rank of JA. The elements in JA will have less additive and multiplicative 
terms than the elements in J and it is thus more suitable for evaluation than 
J. 

However, in JA we still have one inverse matrix left: 2"), Inverse matrices 
tend to have many large sized elements, i.e. many additive and multiplicative 
terms, which make Jacobian matrices containing inverse matrices not very 





well-suited for automatic evaluation. However, if it holds true that the 
restrictions on Q imply the same restrictions on 9", then 27 can be 
reparameterized easily, without this affecting the rank of JA. For example if 
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2Q= : ‘ , where @ and $¢ are symmetrical or diagonal, then 27’can be 
reparameterized like 2: % = 27. A computer program could be written to 
construct the resulting Jacobian matrix, which could be evaluated by ERA or 
Maple. However, even this Jacobian matrix still has many multiplicative and 
additive terms (e.g. A'2"A), and also a large dimension, which makes this 
Jacobian still not very well-suited for automatic evaluation. Furthermore, 
this Jacobian matrix is not generally applicable for all LISREL models, since 
the restrictions on § do not always imply the same restrictions on 27. 

The following theorem (Bekker, personal communication) provides a general 


approach to remove inverse matrices from a Jacobian matrix. 


Theorem 6.4.1 


Let B be a non-singular nxn—matrix. The partitioned matrix 


X, X, Xz 
X = | X, AB'C X,g 
X, X¢ Xq 
X, X, 0 X; 
: ; F ’ X, 0 A Xe 
is of full column rank if and only if the partitioned matrix 
0 -C BO 


is of full column rank. 


Proof 

X%, By 8 
xX, ACO X 
0 0 I, 0 
x Xe ft X 


Define X¥ = , then we find that rank(X ) = rank(X)+n. Now 


consider XZ ZZ, where Z,, Z, and Zz are conformably partitioned matrices 


I00 oO I 0 0 1000 
0 I ABO 0 0 0 0100 
4, = » 242 = , 23 = 
oo! oO 0 CIO 008B 0 
000 TF 00 o7F ooo) 
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We find 
x, ke 2% 
Rag. = X, 0 A X& 
0 -C B 0 


Ge Me Re 


As the matrices Z,, Z, and Z3 are non-singular, rank(XZ,Z2Z3) = rank(X) = 
rank(X) + n. Consequently X is of full rank if and only if XZ,Z2Z, is. Q.e.d. 


We can now apply Theorem 6.4.1 twice to transform (6.4.11) into 


(1,24+Ki,1) (Ae! 2), -(1j2+Ky ,1) (A@A), 0 0 (A@A) I,2 
-~Agl, 0 Nel, 0 0 0 
Fe 0 -A@A 0 Nel, 0 0 
0 0 A’el, 0 0 0 
0 0 R’ 0 A'e@A 0 0 
0 0 0 0 AeA 0 
0 0 0 0 0 T)2 
(6.4.12) 


This Jacobian matrix does not contain any inverse matrices, so that, after 
construction, this matrix can be evaluated by ERA. In principle, a computer 
program could be written that automatically creates this Jacobian matrix, 
given the parameter matrices. We have not developed such a program, because an 
even better approach is possible. This approach has been suggested by Bekker 


(personal communication) and will be discussed in the next section. 


6.5 A relatively simple rank condition 
Let Y, A, A, 2 and = be defined as in (6.4.1) and (6.4.2): 


S=AA'NA'A's E, (6.4.1) 





where Z, 2 and = are the covariance matrices of (y’,x'), (¢',€’) and (€’,6), 


resp., and 
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a= [4 e 


a=[? "|, on |e | (6.4.2) 
0 A 7 


o 1 0 %& 





Let the linear restrictions on the elements of =, A, A and 2 be given by 


vec(=) 
vec( A) 
vec(A) 
vec(2) 


=T, (6.5.1) 


where the symmetry restrictions on = and 92 are not taken into account. 


The symmetry restrictions on 92 are 
2-2 = 0, or (1p2—Kyy) vec 2 = 0. (6.5.2) 


The matrix I,2—K;,,, contains dependent rows: 2 is a kxk-matrix, so we have 
¥ek(k-1) independent rows,  {(e,;@e;)'— (€;@e;)’ | i>j}, where e; is a column 
vector where element i equals one, the remaining elements equal zero. Let the 
matrix S consist of exactly these ¥%k(k-1) rows. So, in addition to (6.5.1), we 


have 
S vec 2 = 0. (6.5.3) 


If 2 is restricted to be symmetric then the symmetry of = need no longer be 
imposed, since it is implied by the symmetry of LY. On the other hand, since 2 
is restricted to be symmetric, the equations (6.4.1) can be written as 


vec © = vec Y — *(I,2+K;,1) vec(AA), (6.5.4) 
where 
A= ANAM. (6.5.5) 


The identifying equations are thus given by (6.5.1), (6.5.3), (6.5.4) and 


(6.5.5), which can also be written as 
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vec( ) —¥2(1,2+K,1) vec(AA) 


R' vec(A ) = 8 
vee) (6.5.6) 
vec(AQA' ) 
S vec(AQA’ ) = 0, 
vec(A-QA’) = 0, 


where $2 can be derived according to Q = AA", 

So (6.4.1) holds and A, A, 2, 2 obey the a priori restrictions imposed on 
the LISREL model, if and only if (6.5.6) is satisfied. To compute the 
corresponding Jacobian matrix, we take the derivatives of the left hand side 
of the identifying equations (6.5.6) with respect to vec'(A), vec'(A), 
vec’(A), vec’(Q), respectively. We get 


—¥e( Ty2tK, 1)(Ti@A) -2(1,2+K; :)(A’@l,) 0 


Re 0 Bei 0 
0 0 ~ Tye 
an 0 0 (1y2+Kyz)(Ix@40) AeA 
0 0 S(142+Ky,~)(Tx@AQ) S( AGA) 
Ty —(1,;92) Kx, 0 -A@I;, 


(6.5.7) 


where we used the fact that, for solutions to (6.5.6), Q=Q'. As 


A’ = AQ, 
(TytKi)(4'@h) = (1e+Ki1)(Ti@A' Ky, and 


Q=A' 2A", 


we may rewrite J as 
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vlTa+Ky (N@A) -¥e(Ty2+Ky,1 )(@AA "2A" ")Ky, 0 0 
R’ 0 I 0 0 
0 0 Iy2 0 
os 0 0 (I, 2+Kp, ¢)x@2A') AGA 
0 0 0 S(A@A) 
de -(11@A"NA'") Ky 0 -A@l, 
(6.5.8) 


When we compute J = VJW with V and W non-singular, and specified as 


I 0 0 0 0 0 
0 IF 0 0 0 0 
0 0 Ff 0 0 0 
V= 0001 (0 0 ; (6.5.9) 
0 0 0 0 Tf 0 
0 0 0 0 0 -Ie@A 
and 
-Iy  -¥AT;@A7*2)Ky, 0 0 
0 %K, (11;@A')Ky, 0 0 
W=l 4g é Red’ 0 , (6.5.10) 
0 0 0 1,@A7 
we get 
A T2+Ky1)(1@A) 0 0 0 
R' 0 %(A’ @!I,) 0 0 
. 0 0 1,@A' 0 
ae 0 0 (Ti2tKe, x) (1,92) Aol, J |f6>1) 
0 0 0 S(A@l;) 
(I, @A) (1;@2) Ky, 1 0 A@l, 
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As V and W are non-singular, we may state that a general LISREL model 
containing linear restrictions on the parameter matrices as in (6.5.2) is 
identified, if and only if J is of full column rank. 

The Jacobian matrix (6.5.11) consists of single elements only: no 
complicated additive and multiplicative terms. This makes it much better 
suited for evaluation than any other matrix discussed in this chapter. The 
only drawback in evaluating this Jacobian matrix is the large dimension: it 
contains (k°+kl +n,) Tows and (2ki+ 2k”) columns. How these dimensions limit the 
evaluation of the identification situation of practical problems is yet to be 
determined. A computer program -— ERALIS - has been written that constructs 
this Jacobian matrix, given the parameter matrices (and their restrictions). 
The rank of this Jacobian matrix can be checked by ERA, or by a CAS, like 
Maple. Example 6.5.1, which is given below, will demonstrate the use of ERALIS 
for a specific model. 

Before discussing the example we shall first, briefly, discuss some 
necessary conditions for identification in LISREL models. These conditions 
provide simple counting rules that have been implemented is ERALIS. These 
rules have been described by Bekker (1990). He gives 24 order conditions for 
the general LISREL model. That is, if one of the conditions given below is not 
satisfied, then almost all points satisfying (6.5.2) will not be locally 
identifiable. The conditions are necessary; they are not sufficient for the 
identifiability of almost all parameter points. 

The conditions are given in tables 6.5.1 and 6.5.2. For each condition 
separately it is assumed that the restrictions on the parameters considered in 
the condition are separable from the remaining restrictions. That is, if the 


restrictions (6.5.2) can be formulated as 


Pe = [® 0 llr | ie (6.5.12) 
0 R be 


then the restrictions of », are said to be separable from the restrictions on 
2. With respect to the conditions in Table 6.5.2 some additional assumptions 
are indicated at the bottom of the table. The conditions give the minimal 
number of restrictions that should be operating on the non—duplicated elements 
of the parameter matrices, or vectors, under consideration. They also give the 
difference between the number of non-duplicated elements in the matrices, or 
vectors, and this minimal number of restrictions. The latter quantity can be 


considered as a maximum for the number of free parameters that can be present 
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in the matrices, or vectors. 

These order conditions are incorporated in ERALIS, as a simple order 
condition check, in order to eliminate at an early stage all non-identified 
models. 


Table 6.5.1 Necessary conditions for identifiability. 


Cond. Parameters Min. # Restrictions Max. # Free Pars 









r(m—p) 










r rows of (B,I,%,) r(mtn-“r+h) 


r rows of (Ay,§%,) r{m—min(n, q)} r{ptmin(n, g)-*r+%} 
r rows of (A,,%5) r{n—min(m, p)} r{qtmin(m, p)-“r+%} 


Br. @ m?+mn—- vam (m+1 ) + 
res min(mn, pn,mq,pq) min(mn,mq, pn, pq) 
Yam (m+1 )+m?+mn— ; 
BI, $;,9, min(mn, mq) ap(p+1)+min(mn,mg) 


vam (m+1)+m?+mn+mp- 
BLP, Ay 8c ,%e min( pn, pq) 






vap(p+1)+min( pn, pq) 






Yon (n+1 ) +m? 4+mn—- ym(m+1)+%q(q+1)+ 
B,D, ©. G6 59s min(mn, pn) min(mn, pn) 















van (nt+1)+m2+mning— |vm(mt1)+49q(q+1)+ 
B,T, A,B, ,Pe,Fs5 min(mq, pq) min(mq, pq) 


BLP, ¢,6_,6,,55 |vm(m+1)+¥an( n+l) +m? 
piace tome aa tag Ln pas 
By tate tt |PmmeL at fapcpo eta ore 


It is assumed that the restrictions on the parameters considered, in each 


condition separately, are separable from the remaining restrictions. 
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Cond. 


It is assumed that the restrictions on the parameters considered, in each 
condition 
Additional assumptions apply to subsets of the Conditions 16-24. 

W.r.t. Conditions 17-19, it is assumed that the restrictions on $, 
separable Condition 
restrictions on the relevant r rows of ®, should be separable. 
W.r.t. Conditions 21-23, it is assumed that the restrictions on 95 
for Condition 20 
restrictions on the relevant r rows of $s; should be separable. 


W.r.t. 


separable from the 


Identification in LISREL models 


Necessary conditions for 


restrictions. 


Parameters Min. # Restrictions 
Ya( r-pt+m)(r—p+m-1)- 


r rows of Ay, 
r{min(n,q)} 


if r> p-m >0 


BT’, Ay,® 2 , 
ot Myo if pom 2m?+mn-min( pn, pq) 


BP, Ay,Oe Pe ,@ 
yore Crt eres 2 = 
eet if pom 2m24+mn-pnt+en(n+1 ) 
BT, Ay hee 3%5,| 2m?4+mn+an(n+1)+ 
ia ai ae eee) ctor Pq 


Perry ers qtn) (r-gtn-1)- 
r{min(m, p)} 


Tr rows Prrow of Ay, Ay, 
r> gq-n >0 


B,',A,,8¢,¢, 


if gon|™ 2in24+mn-min( pn, pq) 


BP, Ag Oc 18e 18 if “% n| vm (m+ )+m?24n?24mn—mqg 


BT, Ay Axe 5% +P 
ae a 


&’ q>n 


» [vm (m1 )+m?24n24mn+ 
is  ——: Pq 


BT BT Ay he BesBp 


2472 3 
f’ psm 2m7+n74mn -pq 





separately, 


from the remaining restrictions; for 


remaining _restrictiotis; 


Condition 24, the restrictions on $ and 95, 


separable from the remaining restrictions. 
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are separable from the remaining 


should both be 


under additional 


Max. # Free Pars . 











r{mtmin(n,q)}- 
2(r—p+m) ( r—p+m-1) 


Yom(m+1 )+m( p—m)+ 
min(pn,pq) 


Yam(m+1 )+m( p-m) + 
Yeq(qt+l )t+pn 


Yam(m+1 )+m( p-m)+ 

fel 
ect | 

¥e(r—qtn) ( r-gqtn—-1) 


von(n+1 )+n(q—n)+ 
ven( m+ 1 )+min(mq , pq) 


van(nt+1 )+n(q—n)+ 
Yep(ptl )+mq 


van(n+1 )+n(q—-n)+ 
eddies 


iene | )+m( p-m)+ 
von(n+1)+n(g-n)+pg 


restrictions. 


are 
16 the 
are 
the 
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Example 6.5.1 
Let a model be specified, like in example 6.3.1, as 


1 Py2 ™m Yii Yi2 éi qa 
+ = 5 (6.3.11) 
Ba, 1 Ne 0 60 £2 2 
1 A, 0 ey 
Yo Az 0 mN &2 
¥Y3 | =] Ag 0 : + | €3 |, (6.3.12) 
Ya 0 A4 Ne E4 
Ys 0 As Es 
xy 1 0 é, 6, 
xq | = [Ag + | @t (6.3.13) 
X3 0 1 §2 5, 


with the variance of 6; = 0, so that 6, = 0 almost surely. Also it is assumed 
that -, &, and ; are diagonal matrices, and ®¢ is unrestricted. 

We will use the program ERALIS to investigate the identification situation 
of this LISREL model. ERALIS needs only as input the parameter matrices and 
the covariance matrices. When given to ERALIS, it yields a sparse 110x96 
Jacobian matrix. We will not give this matrix here, since it yields no 
additional insights. When this Jacobian matrix is evaluated by ERA or Maple, 
we get, just as in example 6.3.1, a null-space of dimension 2. So, as 
expected, the conclusion will be given that the model is not identified. 
However, contrary to the method used in example 6.3.1, we do not get an 
indication of which individual parameters are identified. This is a drawback 
of the method used here. 


If we restrict A, as in example 6.3.1, 
Ay = 4 = 1, (6.5.13) 
then construction by ERALIS yields a 112x96 Jacobian matrix and evaluation by 


ERA or Maple gives no null-space. Consequently, the model, with the additional 


restrictions on A,, is identified. 
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If, instead of the restrictions (6.5.13) we add the restrictions 
Bi2=B ar, (6.5.14) 
and 
A, =1, (6.5.15) 


then construction by ERALIS yields a 112x96 Jacobian matrix. Evaluation by ERA 
or Maple gives once again no null-space. So consequently, this model is 
identified. 

These last two examples can easily be handled by ERALIS. No further manual 
work is needed as was the case in example 6.3.1. However, if the method as 
described in algorithm 6.3.1 can be used, it would be preferable to use that 
method, since information about the identification of the individual 
parameters can be found if the model is not identified. Additionally, the 
evaluation of the Jacobian matrix of this latter method tends to be somewhat 
faster than the evaluation of the Jacobian matrix constructed by ERALIS. On 
the other hand ERALIS has a more general applicability. 


6.6 Conclusion 


In this chapter we have seen a number of practical methods to evaluate the 
identification situation of LISREL models. Two of these methods have been 
implemented using computer algebra. The first method — LISPARAM - is limited 
to evaluating the identification of LISREL models with zero restrictions on 
the elements of B and A, only. The method used in this method is based on 
algorithms described in chapters 4 and 5. The second implemented method to 
evaluate the identification situation of LISREL models — ERALIS - uses a 
Jacobian matrix whose elements consist of single terms only, ie. no additive 
or multiplicative terms, although the number of elements may be large. The 
program can be applied to a variety of models, since the parameter matrices 
are only assumed to be subject to general linear restrictions including 


non-zero restrictions and equality restrictions. 
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CHAPTER 7 


FORMULA MANIPULATION IN STATISTICS ON THE COMPUTER: 
EVALUATING THE EXPECTATION OF HIGHER-DEGREE 
FUNCTIONS OF NORMALLY DISTRIBUTED MATRICES 


7.1 Introduction 


In this chapter, based on a paper of Merckens and Wansbeek (1988), we consider 
the case of the evaluation of the expectation of higher—degree functions of 
normally distributed random matrices. We first indicate the relevant 
statistical context, and next show how computer algebra comes in. 

If x ~ N(0,1), it is well-known that 


E(x") = (n-1)(n—3)...3.1 = nt/{2™?(n/2)!}, (7.1.1) 


for n even, and zero for n odd. When the expectation involves a higher—degree 
function of a vector or matrix variate (rather than a scalar), the expressions 
become quite a bit more complicated, as is apparent from the pertaining 
literature. 

For example, if x is a vector of k normally distributed elements with zero 
expectation and covariance matrix V (kxk), so x~N,(0,V), Magnus (1978) 
discusses E(x'Ax)" and, more generally, E(Il7.,x’A;x), where Aj,.--,4y are 
non-random kxk-matrices, and he introduces the complicated concept of 
A(s)-polynomials to characterize these expectations. Further elaborations were 
subsequently given by Don (1979) and Magnus (1979). 

As to normally distributed random matrices, Neudecker and Wansbeek (1983), 
for example, contains a number of results for fourth—degree functions, like 
var(vec(X@X)); X is a matrix with var(vecX) =V, with V some non-random positive 
semi-definite (p.s.d.) matrix. They also consider E(X@X@X@X) for the more 
restricted case where V=WeI, with W some non-random p.s.d. matrix and J the 
identity matrix. More results are contained in Neudecker and Wansbeek (1987) 
(hereafter denoted by NW). E.g., if E(X)=M and var(vecX)=VeU, U some 
non-random p.s.d. matrix, NW consider the covariance between vecX’AX and 
vecX'BX, and E(X’AXCX’BX), with A, B, C non-random matrices, not necessarily 
p.s.d. or even symmetric. To get an idea of how complicated things become, 


consider the latter: 
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E(X’AXCX'BX) = {(trAUV+M'AM}C{(trBUV +M’BM} + 
(trAUBUVCV + VC’M'A'UBM + M'AUB'MCV 4 
(trAUBU)(trCV V + (trAMCM'BUV + (trCV)M’AUBM, (7.1.2) 


where "tr” stands for the trace of a matrix. To prove this result, NW develop 
and use a method that we might call "repeated conditioning”, and that works as 
follows. Label the four X-matrices from 1 through 4, and rewrite the 


E-operator as 
E=E,2E3q + Eyskoq + ExaEss, (7.1.3) 


where E,;, t,j=1,...,4 is meant to denote the expectation with respect to X; 
and X;, considering the other X’s as constants. It is useful at this place to 
consider the relationship between (7.1.2) and (7.1.3). For didactical reasons 
it suffices to look at the simple case U=I, V=I, M=0. We use the results 
E(X’PX) = (trP)I, E(XQX) = Q', for non-random matrices P, Q of appropriate 
order. Then 


E(X'AXCX'BX) = 
= (Ey2E3q + EygEgq + Ey4E3)(Xy'AX2CX3'BX4)= 
= Eyo(Xy'AX,C)(trB) + Eyg(Xy'AB'X3C’) + E,4(X{'ABX,)(trC) = 
= (trA)(trB)C + (trAB’)C’ + (trAB)(trC)I, (7.1.4) 


and this is exactly (7.1.2) for this simple case. This method of evaluating 
the expectation of fourth-degree functions of random matrices is quite 
straightforward, though restricted to the normal case. 

The purpose of this chapter is twofold. First, it is shown in section 7.2 
that a recent result due to Guiard (1986), makes it possible to extend the 
NW-method to functions of degree higher than four. Second, it is argued in 
section 7.3 that, for any given problem instance, the analytical evaluation of 
expectations like (7.1.2) is eminently suited for being performed by computer. 
How to do this is also shown. In section 7.4 we will discuss the concept of 
A(s)-polynomials as given by Magnus (1978) and Don (1979). We will show that. 
applying our method for a specific subset also yields their findings. Examples 
of “automatically” generated formulas are given in the Example 7.3.1 and the 
Appendix. These formulas can next be transformed into a source text in a 
higher-level programming language like Pascal. An example of such a text is 


also given. 
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7.2 Evaluation of expectations 


Consider the set T, ={1,2,...,n}, m even. Using each element once, we can 
construct a set containing unordered pairs of elements from T,, e.g. 
{(1,2),(3,4),---,(m—1,n)}. This set contains n/2 elements, and many such sets 
can be constructed: there are clearly (n—1).(n—3)...3.1 different ones. This 
is precisely the right-hand side of (7.1.1), and it provides a justification 
for a breakdown a la (7.1.3) when we have to do with a standard normal scalar. 
As another example, for n=6 there are (6—1)(6—3)(6—5) =15 such sets, and the 
corresponding breakdown of the expectations operator then reads as 


E = EyQEggEsg + ExcEssEge + Exok sels 
Ey3Eo4Ese + EqsEosEyg + Ey3E26k4s 
EysEggEsg + EysEosEse + ExaE ook ss 
EysEo3E46 + EjsEogEs6 + EisE26E34 
EygEosEqs + ExgEoaEss + EyeEasE 3a. 


(7.2.1) 


++ 4+ + 


So the repeated conditioning method gives a means to evaluate E(x”) with n 
even and x standard normal. As the outcome is known, this method is useless in 
this restricted context, but would be very useful if it were applicable in 
more general situations, in particular the following one. Let %j,...,x, be 
jointly normally distributed with mean zero and variances and covariances 9;,;, 


we would like to evaluate 
n 
E( 1 x;). (7.2.2) 


That repeated conditioning also applies here follows from a result due to 
Guiard (1986). Define the function M(R,) recursively as 


M(Ro) = 1, 
M(Ry) = Lier,,\(e)7HM(Rp\{k,1}), Pp = 2,4,...,n. (7.2.3) 


where R, denotes any subset of T,, containing p elements (cf. Guiard, 1986; pp. 
280-281). The choice of the index k is arbitrary. 
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Theorem 7.2.1 


em Ve the moment generating function of 


Let ¢ be an m-vector and W(t) = 
the multivariate normal distribution with expectation zero and covariance 


matrix V = {o,;}. Let |x] be the entier function of %y H (Xje1 ta) = 1, 
ERy 
and M(R,) defined as in (7.2.3), then 


a w(t) 
Bt, dtg,-.., Ok, 


Ln/2] . 
W(t) Dae bale (Sh ae | M(T,-R | (7.2.4) 


n-2q 


For a proof of Theorem 7.2.1 see Guiard (1986; pp. 284-287). 


Theorem 7.2.2 
Let x,,...;X, be jointly normally distributed with mean zero and variances 


and covariances 0;;, and M(T,) defined as in (7.2.3), then 


E( fx.) = M(T,), n even, (7.2.5) 
E( fh x4) = 0, n odd. (7.2.6) 
Proof 


Using the moment generating function we get 


i _ 9 Vt) x 
E(*i) = 35, dtg,.0;0te| (t=0) = 


Ln/2] 
v(0) | ae (Zia1 tu) MCTy-Re 


q=0 n-2q (t=0) 


Since the term TH (Xj. toy) # 0 only for n-2g = 0, this reduces to 
TeRn_2g 


E( I x,;) = b, ~ M(T,-R,) = M(T,), for n even, and 


On 


E( fix) = 0, for n odd. 


i= 


= 7 114 











Chapter 7 


For example take the simple case n=4 (and k=1 in the first step), then 
evaluating (7.2.3), yields 
M(Ry) 012M(3,4) + O13M(2,4) + 644M(2,3) = 


12034 + 13% 24 + %14%23, (7.2.7) 


which yields the same result as from replacing E in (7.2.2) by (7.1.3). 

So, for the scalar case, repeated conditioning is a valid method, and it 
directly extends to the case where the expectation of a matrix is involved. 
Consider for example E(Y) with Y = X'AXCX'BX, with vecX normal with zero mean 
and arbitrary covariance matrix. This is easily reduced to a number of scalar 


problems by using the equivalence: 
trE(YP) = trQP for any P + E(Y) = Q. 


Now trE(YP) is just a sum of terms like (7.2.2), with n=4 and xj,...,%4 for 
each term denoting some element from one of the four X’s in X'AXCX'BX each. 
Due to the linearity of both summation and expectation, we might as well apply 
(7.1.3) to the problem in matrix format directly. This was done by NW to 
arrive at (7.1.2), but there a specific argument for n=4 was used. The Guiard 
result shows that matrix problems for any even n can be handled in the same 
way. 

The potential for the method stretches beyond problems like E(X'AXCX'BX) 
and its direct extensions. First, note that the pattern of transposed and 
non-transposed X’s can be arbitrary, second, if X has a non-zero mean, 
adaptation is straightforward, for n even and for n odd. Third, NW show that 
cov(X'AX,X’BX) can be directly obtained by manipulating on E(X’AXCX'BX), so 
repeated conditioning is also useful here. Evaluating this covariance has been 
considered by other authors by a variety of methods, e.g. based on the Wishart 
distribution, following the results by Magnus and Neudecker (1979), by 
Neudecker (1985); on moment-—generating functions - Von Rosen (1987); or on 
matrix differentiation methods applied to cumulant generating functions -— 
Browne and Neudecker (1988). But the method of repeated conditioning is simple 


and general, and extends to higher-degree moments straightforwardly. 
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7.3 Symbol manipulation by computer 


From now on, we consider the ”generic” problem of evaluating 


n/2 
ef II X'P,XOQ; ) (7.3.1) 
t=1 





where n is a given even integer, and the P,’s and Q,’s are non-random matrices 
(which may or may not be symmetric); X is normally distributed with E(X) =0 and 
var(vecX) = V @U. The two distinctive features of this problem are that (a) 
the resulting expression is long and complicated, and (b) going through the 
derivation by means of repeated conditioning is a tedious affair but takes 
place by repeated application of a small number of "production rules”. Hence, 
the job is eminently suited for being performed by computer. We have developed 
such a program, which runs on a PC. The construction of a program like this — 
ie., able to manipulate analytical formulas — can be handled in several ways. 

The easiest way is to use a symbol manipulation language like Reduce or 
Macsyma. By using one of these languages, problems like ‘how to recognize and 
store formulas’ or ‘how to multiply two formulas’ have been taken care of. 
Sometimes, however, these languages do not provide enough possibilities to 
tackle the problem, or these languages are not at hand. Then the basic 
problems like the ones mentioned have to be dealt with by, for example, using 
Pascal and by defining tree structures for storing and manipulating formulas. 
In chapter 2 some basic information is given how one could proceed. However, 
the examples given in that chapter are not straightforwardly transformed in a 
computer algebra program capable of handling (7.3.1). A new abstract data type 
(see section 2.4) has to be defined: a data type capable of holding the 
necessary information on matrices and traces, and algorithms to manipulate the 
information. Since available computer algebra systems are not able to handle 
our type of matrix problems directly, we opted for Pascal. In this chapter we 
will not discuss the technical aspects of the implementation of this computer 
program, we only show which algorithms are needed. 

Our program consists of four parts. The first part generates all 
*two-partitions” like (7.1.3) or (7.2.1); the second part contains the 
"production rules” needed for evaluation of expectations over two X’s; the 
third part recognizes and collects identical terms thus spawned and in the 
fourth part the result obtained is translated into a direct-to-use computer 


program. 
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7.3.1 Generating the two-partitions 


Given dimension n, the number of partitions is simply generated by the 
recursive procedure (7.2.3). Since in the programming language Pascal one can 


use sets, this formula can easily be programmed. 


7.3.2 The production rules 


When for the given n the set of two-partitions is derived, it remains to apply 
the E-operator, broken down into the (n-1)(n-3)...3.1 additive terms each 
consisting of n/2 expectations involving two X’s. Each of these expectations 


has one of the following forms, which are presented together with their 


outcome: 
E(trX’AXB) = trAU.trBV (7.3.2) 
E(trXCXD) = trUCVD (7.3.3) 
E(trFX)(trGX) = trF'VGU (7.3.4) 


These rules (for proofs see NW), together with the analogous results where all 
X’s are transposed, constitute the production rules. Note that in any 
particular step in some application, the matrices A, B,.C, D, F, G are in 
general functions of the P,’s, Q;’s, U, V and of those X’s that are considered 
constant at that moment. ; 


7.3.3 Collecting identical terms 


When all P,’s and Q,’s are different and non-symmetric, all (n—1)(n-3)...3.1 
terms from (7.3.1) will be different. As soon as some kind of structure is 
present, some terms become identical, and should be collected. 

To allows the program to recognize identical terms, all terms have to be 


ordered in a unique way, so rearrangements like 

tr(ED)tr(CB’A’) = tr(ABC’)tr(DE) (7.3.5) 
have to be performed. As can been seen from this simple example, two kinds of 
ordering have to be defined, (a) an ordering within a particular trace (the 
"trace ordering”), and (b) an ordering of the traces within a particular term 


(the "factor ordering”). 
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(a) The trace ordering 


Two basic rules concerning traces are tr(AB) = tr(BA), and tr(C) = tr(C’). So 
when the trace is taken over the product of k (say) matrices, there are 2k 
ways to represent it. The ordering that we use is the lexicographic ordering 
based on the ASCII table, with the understanding that symbols with 
transpose-indicator come in the ordering immediately after their 
™untransposed” counterparts. For example, tr(P,'Q2'P;P2:'Qy'), tr(P1Q,P2P1'Q2) 
and tr(Q,P,P,'Q,P;) will all be standardized” to tr(P,P2'Q,'P)'Q,’). 


(b) The factor ordering 


Since of course tr(A)tr(B) = tr(B)tr(A), an easy sorting to establish a unique 
ordering is as follows. Let #A stand for the number of matrix factors in A. 
Then if #A>#B, put tr(A) in front; if #B>#A, put tr(B) in front; if #A=#B, use 
the lexicographic ordering. 

Now a unique ordering has been reached and the computer can easily compare 
terms and collect identical ones. A few examples are given below; more 
examples can be found in the Appendix. 


Example 7.3.1 
The evaluation of E(tr(X’P,XQ,X'P,XQ.X'P3XQ3)), with E(X) = 0, gives 


tr(P,UP,UP,U) tr(Q,V) tr(Q.V) tr(Q3V) 
tr(O,VOVOV) tr(PU) tr(PU) tr(P3U) 
tr(P,UPU) tr(QVOV) tr(P3U) tr(QV) 
tr(P,UPU) tr(QVOV) tr(PV) tr(Q,V) 
tr(P,UPU) tr(QVQV) tr(P,U) tr(QV) 
tr(P,UP,UP3'U) tr(QVQ3,V) tr(QV) 
tr(P\UP'UP3U) tr(QVO.V) tr(Q,V) 
tr(P,UP;'UP,'V) tr(QVO,V) tr(QV) 
tr(Q,VOV0,V) tr(P,UP;'U) tr(PU) 
tr(OVO,VO,V) tr(P2UP3'V) tr(P,U) 
tr(Q,VO,'VO,V) tr(P,UP,'U) tr(P3U) 
tr(P\UP,'UP3'U) tr(QVQVQ.V) + 
tr(P,UP,UPU) tr(QVOVOV) + 
tr( 
tr( 


ttt tt 


+etryett 


P,UP,UP,U) tr(Q,VO;VOV) + 
P,UP;UPV) tr(QVQ,VQ3V). 
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Example 7.3.2 
The evaluation of E(tr(X’P,XQ,X'P,XQ,X'P,XQ3)), compared with examp-‘e 
7.3.1: Pj=P,=P3, Q,=Q2, E(X)=0, gives 


tr(P,UP,UP,U) tr(Q,V) tr(QV) tr(QV) + 
tr(Q,VOVOsV) tr(PwW) tr(P,U) tr(P\U) + 
tr(P,UPWU) tr(Q,VO,V) tr(PW) tr(Q,V) + 
2 * tr(PUPU) tr(QVOsV) tr(PU) tr(Q,V) + 
tr(P,UP,UP,'U) tr(O,VO,V) triQV) + 

2 * tr(P\UP,UP,'U) tr(OVO,V) triOV) + 
tr(O,VO,VO,'V) tr(P\UP,'U) tr(PU) + 
tr(Q,VO,'VOnV) tr(P\UP,U) tr(PU) + 
tr(Q,VOVO,V) tr(PyUP,U) tr(PW) + 
tr(P\UP,UPU) tr(QVOVOV) + 
tr(P,UP,UP,U) tr(QVOVO,V) + 
tr(P,UP,UP,'U) tr(Q,VO;VOV) + 
tr(P,UP,UP,'U) tr(Q,VO,VO,'). 


Example 7.3.3 
The evaluation of E(tr(X’P,XQ,X'P,XQ,X'P,XQ,)), compared with example 7.3.1 
P,=P,=P3, Q;=0,=Q3, E(X)=0, P; and Q; symmetric, yields 


tr(P,UP,UPU) tr(Q,V) tr(Q,V) tr(QV) + 
tr(Q,VOVOV) tr(PU) tr(PU) tr(P\U) + 

3 * tr(P,UPU) tr(QVOV) tr(PU) tr(QV) + 
3 * tr(P,UP,UP,U) tr(QVOV) tr(QV) + 

3 * tr(QVOVOWV) tr(P\UPU) tr(PU) + 

4 * tr(P\UP,UPU) tr(QVOVO,V). 


Example 7.3.4 
The evaluation of E(tr(XP,X’Q,X'P,X'Q,XP3XQ3)), compared with example 
7.3.1: the X-matrices are differently transposed, yields 


tr(PVP3U) tr(QVO3U) tr(P,V) tr(O.U) + 
tr(P\VO,'UPVP,U0,V) tr(Q,U)+ 
tr(P,VO,UPVPUQV) tr(Q,U)+ 
tr(P.VO,'U03'VP,UQ,'U) tr(P,V)+ 
tr(PVO,UO,VP,UOU) tr(P,V)+ 
tr(P\VP\UOUOV) tr(PVOU)+ 
tr(PVQ,UO.UP,V) tr(P3U0V)+ 
tr(P,VP,'UP3V) tr(Q,VO3U0,'U) + 
tr(P,VP3UPV) tr(QVOU0U) + 
tr(P,VP,U0Q,V) tr(P,VO,'UO,'U) + 
tr(P,VO,'UP,V) tr(P3U0,'U0,V) + 
tr(P,VP,'U0,'U0,VP,UO,V) + 
tr(P,VP_UQ;'VP,UQ,'UOV) + 
tr(P,VQ,'UPVO,U0,'UP3'V) + 
tr(PVOU0,'UPVO,'UP;). 
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7.3.4 Generating a computer program 


As an integral part we have included an option to translate the ‘freshly’ 
computed formula into a Pascal computer program that is able to actually 
compute the desired expectation of given numerical values of the various 
matrices. So the tedious job of programming cumbersome formulas like those 
presented in the Appendix can also be taken over by the computer. Example 
7.4.4 in the Appendix shows the result.of transforming example 7.3.3. 

Various extensions are conceivable. One is to make the pattern of 
transposed and non-transposed X’s arbitrary. This is a straightforward 
generalization and can be implemented without much effort; see example 7.3.4 
to give an idea. Another extension that has been implemented is to relax E(X) 
= 0 and have E(X) = M instead. By substituting X = Y+M, E(Y) = 0 and 
elaborating the product, 2” new problems arise that all involve 
zero—expectation matrices as before; see example 7.A.1. Example 7.A.2 shows 
that the number of X-matrices has not necessarily to be an even number; of 
course all terms then vanish when M=0. Also, the structure U@V on the variance 
of vec(X’) can be relaxed, and can be replaced by- an arbitrary symmetric 
matrix instead. This requires a generalization of the production rules (7.3.2) 
- (7.3.4). This extension has not been implemented at present. The current 


k nj 
implementation can handle all problems of the form B(er( I ZiyPsy)), with 
Zjz=X or Z;j=X', E(X) = M or 0, and any P,; can be symmetric, equal to P,, or 
equal to an identity-matrix. The information given will be used to simplify 


the input expression and the resulting expectation. 


7.4 Overlap with A(s)-polynomials 


As mentioned in the introduction, Magnus (1978, 1979), and Don (1979) 
discussed in three papers E(Il;.,x'Ajx), where x~N,(0,V) and Aj,...,An, are 
non-random symmetric kxk—-matrices. Note that the symmetry assumption on A, is 
non-restrictive: if A; is non-symmetric, we replace x'Ajx by x'A;x, with A; = 
A,A; / 2. To characterize these expectations Magnus introduces the concept of 
A(s)—-polynomials. Don formulates a general formula for the A(s)—polynomials. 
We will suggest in this section that applying our method to (7.3.1) for 


particular choices for some of the matrices yields exactly an A(s)—polynomial. 
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Chapter 7 


7.4.1 The A(s)-polynomial 


Let Aj,A,,...,4, be real symmetric matrices of the same order. The following 
definitions, given by Magnus, describe the A(s)—polynomial. 


Definition 7.4.1 A(s)-form 

Divide the index set {1,2,...,s} into mutually exclusive and exhaustive 
subsets. Within each subset take the trace of the matrix product of the A,’s 
corresponding with indices from the subset. The product of all these traces 


will be called an A(s)—form. 


Example 7.4.1 
tr(A,)tr(A,A3),tr(A;)tr(Az)tr(A3), and tr(A,A,A43) are all examples. of an 
A(3)—form. 


Definition 7.4.2 Similarity class 
Two A(s)-forms belong to the same similarity class if their corresponding 
subsets differ only by a permutation of indices. 


Example 7.4.2 

tr(A,)tr(A,A3) is equal to tr(A3A,)tr(A,), but not necessarily equal to 
tr(A)tr(A,A3), However, all these A(3)-forms belong to the same similarity 
class. On the other hand tr(A,A,A3) belongs to a different similarity class. 


Definition 7.4.3 A(s)—sum 
The sum of all unequal A(s)—forms within a similarity class is called an 


A(s)-sum. 


Example 7.4.3 
The three A(3)-sums are: tr(A,)tr(A,43) + tr(A,)tr(A,A3) + tr(A3)tr(A,A), 
tr(A,A,A3), and tr(A,)tr(A2)tr(A3). 


Definition 7.4.4 A(s)- polynomial 


Any linear combination of A(s)-sums is called an A(s)—polynomial. 
Example 7.4.4 
The A(3)-polynomial is v,(tr(A,)tr(A,A3)+tr(A2)tr(A,A3)+tr(A3)tr(A,A2)) + 


vgtr(A,AgA3) + vgtr(A,)tr(A,)tr(A3), with v; constants. 


121 








The expectation of higher—degree functions of normally distributed matrices 


Theorem 7.4.1 (for a proof see Magnus (1978) or Don (1979)) 

Let A;, Ao, ..-, Ag be real symmetric nxn matrices, x ~ N(0,V), V positive 
definite. Then E(Ij_,x’A;x) is an A(s)-polynomial with each A; replaced by 
AW. 

iV. 


Don (1979) shows that the coefficients of an A(s)-polynomial can be 
computed easily by using all possible partitions of the number s. He defines a 
partition A of s as a sequence of the integers A(1) 2 A(2) > ... 2 A(r) 2 1, 


with TA(i)=s, r arbitrary. 


Example 7.4.5 
If s=4 then all possible partitions are (1,1,1,1), (1,1,2), (1,3), (2,2), 
and (4). 


Now for any partition of s define n; (j=1,...,5) as the number of times that j 


appears in the partition. Of course it holds true that Ljn,=s. 


Example 7.4.6 : 
For s=4 the following table holds the information on the partitions and the 


values of n;. 


partition Ny Ng Nz =~ 
1111 4 

112 2 1 

13 1 1 

22 2 

4 1 


Don (1979) shows that the coefficients of the A(s)-polynomials are equal to 
28-"-"2 and that the number of non-equal A(s)—forms in a similarity class 
equals s!2%1+"2 / T3217 j!(2j)"4. 

We can use this result to test if our algorithm is correctly implemented. 
In (7.3.1) we consider E} tr m/? X'P;XQ;|, with n even and X a matrix, E(X)=0 
and var(vecX’) = U@V, whereas Magnus and Don consider E(tr Ij_,x'P;x), with x 
a vector, P; symmetric and var(x’)=V. So, if we choose n=2s, P; symmetric, 
O=In, V=Im, M=1, then these two formulas coincide. So computation of (7.3.1) 
with these restrictions should yield an A(s)-polynomial. Indeed the following 
example shows for n=8 that the result is an A(4)—polynomial. 
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Example 7.4.7 
The evaluation of E(tr(X’P,XX'P,XX'P,XX'P,X)), with E(X)=0, (Q; = J). P, 
symmetric, yields 


ct 


r(P,) tr(P,) tr(P3) tr(P4) + 
tr(P,P.) tr(P3) tr(P4) + 
tr(P,P3) tr(P,) tr(P4) + 
tr(P,P,) tr(P,) tr(P3) + 
tr(P2P3) tr(P;) tr(P4) + 
tr(P2P4) tr(P,) tr(P3) + 
tr(P3P4) tr(P,) tr(P2) + 
tr(P\P2P3) tr(P4) + 
tr(P,P2P,4) tr(P3) + 
tr(P,P3P,4) tr(P2) + 
tr(P2P3P4) tr(P;) + 
tr(P,P2) tr(P3P4) + 
tr(P,P;) tr(P2P4) + 
tr(P,P,4) tr(P2P3) + 

16 * tr(P,\P2P3P4) + 

16 * tr(P,P,P4P3) + 

16 * tr(P,P3P,P4). 


mh RP OOO CON NNNNN 
* %¥ 4% HH KKH HK HH *K 


We have shown in this section that the computation of A(s)-polynomials can be 
done by using our general method for determining the expectation of 
higher-degree functions of normally distributed matrices, for some specific 


choices of matrices. 
7.5 Conclusion 


We have shown in this chapter that computer algebra can be used in the context 
of higher-degree functions of normally distributed matrices. This context may 
appear to be of limited practical importance, but there have nevertheless been 
a fair number of papers in this field, and with the results presented here we 
are able to have the computer reproduce numerous formulas that have appeared 
in the literature as interesting results per se, plus a wide range of 
generalizations. 

In the next chapter we show a computer algebra program capable of deriving 
estimators in complicated ANOVA models. These problems have in common with 
this problem that the derivations in many instances are both highly 
complicated and tiresome to perform by hand, and rather straightforward in the 
repeated application of a limited set of production rules. For this type of 


problems computer algebra can also be a useful tool. 
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7.A Appendix 
Example 7.A.1  E(tr(X'P,XQ,X'P,XQ,X'P,XQ,)) ,P, and Q,; symmetric, E(X)=M 


tr(P,UP,UP,U) tr(Q,V) tr(Q,V) tr(Q,V) + 
tr(Q,VOVOV) tr(P,U) tr(P,U) tr(PU) + 
3 * tr(P,UPU) tr(Q,VO,V) tr(PU) tr(Q,V) + 
tr(MQ,M’P,UP,UP,) tr(Q,V) tr(Q.V) + 
tr(MO,VO,M'P,UP,) tr(PyW) tr(QV) + 
tr(MO,VO.VO,M'P,) tr(P,U) tr(P\U) + 
tr(MO,M'P,UP,) tr(Q,VO,V) tr(PU) + 
tr(MO,VO,M'P,) tr(P,UP,U) tr(QV) + 
tr(P,UP,UP,U) tr(QVQ,V) tr(Q.V) + 
tr(QVOVOV) tr(PUPWV) tr(PU) + 
tr(MQ,M'P,MQ,M'P,UP,) tr(Q,V) + 
tr(MQ,M’P,MO,VO,M'P,) tr(P\U) + 
tr(MO,VO,M'P,UP,UP,) tr(Q,V) + 
tr(MO,VO,VO,M'P,UP,) tr(P\U) + 
tr(MO,M'P,UP,UP,) tr(Q,VO,V) + 
tr(MO,VOVO,M'P,) tr(PUP WU) + 
tr(MQ,M'P,UP,) tr(MO,VO,M'P,) + 
tr(P,UP,UP,U) tr(Q,VO,VO,V) + 
(MQ,M'P,MQ,M'P,MQ,M'P,) + 

6 * tr(MQ,M'P,MO,VO,M'P,UP,) + 

12 * tr(MO,VO,VO,M'P,UP,UP,). 
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If M will be set to 0 then the outcome will reduce to the outcome of 
example 7.3.3 (as it should). 


Example 7.A.2  E(tr(X’P,XQ,X'P2XQ2X'P3)) , E(X) = M (dimension is 5) 


tr(MP,'UP,/UP,') tr(Q,V) tr(Q,V) 
tr(MQ,VP,'UP,') tr(P,U) tr(Q,V) 
tr(MQ,VO,'VP,') tr(PyU) tr(P2V) 
tr(MP3'UP,') tr(Q,VO,V) tr(PU) 
tr(MQ,'VP,') tr(P\UPV) tr(O,V) 
tr(MPs'MQ,'M'P,'UP,’) tr(Q,V) + 
tr(MP3'MQ,'VO,'M'P,’) tr(PU) + 
tr(MP,'UP,'MO,'M'P,') tr(QV) + 
tr(MQ,VP,'UP,UP;) tr(Q.V) + 
tr(MQ,VO,VP,'UP,) tr(PU) + 
tr(MQ,'VP3'MQ,'M'P,') tr(PU) + 
tr(MOVP,'UP,UP,) tr(QV) + 
tr(MO.VO,VP;'UP,) tr(PU) + 
tr(MP;'UP.UP,') tr(Q,V0,V) + 
tr(MQ,M'P,UP,) tr(MQ,VP;’) + 
tr(MQ,VQ,'M'P,') tr(MP3;'UP;') + 
tr(MQ,'VO,VP,') tr(P,\UPU) + 
tr(MP3'MQ,'M'P,'MQ,'M'P,') + 
tr(MP3'MO,VO,M'P,UP,’) + 
tr(MP;'UP,MO,VO,'M'P,’) + 
tr(MQ,M’P,MO,VP,'UP,) + 
tr(MQ,VP3'MO,'M'P,'UP,) + 
tr(MOQ,VO,'VP,UPUP,) + 
tr(MOQ,'VO.VP3;UP,UP,') + 
tr(MQ,VQ,VP,'UP,'UP,') + 
tr(MO,VO,VP,'UP, UP,). 


tte et 
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Example 7.A.3 dimension 10, E(X) = 0, P;=P;, Q;=Q;, P; and Q; symmetric 


E(tr(X’P,XO,X'P,XQ,X'P,XQ,X'P,XQ,X'P,X Q;)) 


tr(P,UP,UP,UP,UPW) tr(Q,V) tr(Q,V) tr(Q,V) tr(Q,V) tr(Q.V) + 
tr(QVOVOVOVOWV) tr(PU) tr(PU) tr(PU) tr(PwW) tr(P\U) + 
5 * tr(PUPUPUPWV) tr(QVOV) tr(PW) tr(Q,V) tr(Q,V) tr(Q,V) 
* tr(OVOVO,VO,V) tr(P,UP,U) tr(PU) tr(P,U) tr(PU) tr(Q,V) 
* tr(P\UP\UPU) tr(QVO,VO,V) tr(P) tr(PU) tr(Q,V) tr(O,V) 
* tr(P\UP\UPU) tr(P,UP,U) tr(QVOV) tr(Q,V) tr(Q,V) tr(O,V) 
* tr(P\UPUPW) tr(Q,VOV) tr(QVO,V) tr(PW) tr(PW) tr(Q,V) 
* tr(OVOVO,V) tr(PUPU) tr(P,UP,U) tr(P,U) tr(Q,V) tr(Q,V) 
* tr(QVO,VO,V) tr(PUPW) tr(QVO,V) te(PV) tr(PW) tr(PU) 
* tr(P\UPU) tr(P,UPU) tr(QVO,V) tr(Q,VO,V) tr(PWU) ix(,V) 
10 * tr(P,UPUPUPUPV) tr(QVO,V) tr(QV) tr(Q,V) tr(Q,V) 

10 * tr(0,VO,VO,VO,VO,V) tr(P,UP,U) tr(PU) tr(P) tr(PU) 
15 * tr(P\UP,UP,UPU) tr(QVOVOV) tr(PW) tr(Q,V) tr(Q,V) 
15 * tr(QVOVO,VO,V) tr(P,UP\UP,U) tr(PU) tr(P) tr(Q,V) 
15 * tr(P\UP,UP,UPU) tr(Q,VO,V) tr(Q,VO,V) tr(PW) tr(0,V) 
15 * tr(QVOWVO,VO,V) tr(P,UPU) tr(P,UP,U) tr(P,U) tr(O,V) 
15 * tr(P,UP,UPU) tr(QVOVO,V) tr(PUPW) tr(Q,V) tr(Q,V) 
( 
( 


aoa cr orci a 
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15 * tr(P\UP,UPU) tr(Q,VO.VOV) tr(QVO,V) tr(P,U) tr(P\U) 
10 * tr(PUP,UPWV) tr(P\UP,U) tr(Q,VO,V) tr(Q,VOV) tr(O,V) 
10 * tr(QVO,VO,V) tr(P\UPU) tr(PUP,U) tr(Q,VO,V) tr(P,U) 
40 * tr(P,UPUPUPUPU) tr(QVOVO,V) tr(Q,V) tr(Q,V) 

40 * tr(O,VO,VOVOVO,V) tr(P,UP,UPWU) tr(PWU) tr(PU) 
25 * tr(P\UP,UP\UP,UP,U) tr(Q,VO,V) tr(Q,VO,V) tr(Q,V) 
25 * tr(OVOVOVOVOWV) tr(PUPU) tr(P\UPU) tr(P,U) 
60 * tr(P\UP\UP,UP,U) tr(Q,VO,VOVOV) tr(PU) tr(Q,V) 
45 * tr(PUPUPUPW) tr(QVOVO,V) tr(OVO,V) tr(PU) 
45 * tr(QVO,VOVOV) tr(PUP,UPU) tr(P,UPWU) tr(Q,V) 
25 * tr(PUP\UPU) tr(QVO,VO,V) tr(P,UPWU) tr(0,VO,V) 
100 * tr(P,UPUP,UP,UPU) tr(O,VOVO,VO,V) tr(Q,V) + 
100 * tr(Q,VOVOVOVOV) tr(P,UP\UP,UPU) tr(PU) + 
60 * tr(P\UP\UP,UP,UPU) tr(0,VO,VO,V) tr(Q,VO,V) + 
60 * tr(Q,VOWVOVOVOV) tr(PUP PU) tr(P\UPU) + 
148 * tr(P,UP,UP,UPUPU) tr(O,VO,VO,VO,VO,V). 


++ttttte¢est 


tet t+t+ts 
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Example 7.A.4 Generated computer program using example 7.3.3 


PROGRAM DIMEN6; 


CONST 
DIMP1 = 3; DIMU = 3; DIMQ1 = 3; DIMV = 3; n = 6; 
Maxdim = 3; {Maximum of all dimensions} 








TYPE 
names = (P1, U, Qi, V); 
vector = array{1 .. maxdim] of real; 
matrix = array[1 .. maxdim] of vector; 
namevec = array[(l .. n] of names; 
boolvec = array[1 .. n] of boolean; 
intvec = array{Pl .. V] of integer; 
BigMatrix= array[{P1 .. V] of matrix; 
VAR 
fact: array[1 .. 6] of real; 
MAT: BigMatrix; 
pos: namevec; {Contains information on position; } 
{eg. for tr(P1 U), pos[1] = P1; pos[2] := U } 
t: boolvec; {Contains information on transpose-indicators} 


totaltrace: real; 
datafile: text; 
dim: intvec; {Contains all dimensions } 


{The following two procedures are not listed since it is identical for every 
computed formula} . 

function computetrace(.... ): real; 

procedure readdata(....); { Reads datamatrix from file } 


BEGIN 
totaltrace := 0; assign(datafile, DATA’); reset(datafile); 
dim[P1] := DIMP1; dim[U] := DIMU; dim[Q1] := DIMQ1; dim[V] := DIMV; 
readdata(MAT[P1], dimP1, datafile); readdata(MAT[U], dimU, datafile); 
readdata(MAT[Q1], dimQ1, datafile); readdata(MAT[V], dimV, datafile); 
pos[1] := Pl; pos[2] := U; pos[3} := Pl; pos[4] := U; pos[5] := Pl; 
pos[6] := U; 
t{1] := FALSE; t[2] := FALSE; t[3] := FALSE; t[4] := FALSE; 
t[5] := FALSE; t[{6] := FALSE; 
fact(1] := computetrace(MAT , dim, 6, pos, t); {fact{1]=tr(P,UP,UPU)} 
pos[1] := Ql; pos[2] := V; t[1] := FALSE; [2] := FALSE; 
fact[2] := computetrace(MAT , dim, 2, pos, t); 
pos[1] := Q1; pos[2] := V; pos[3] := Ql; pos[4] := V; pos[5] := Ql; 
pos[6] := V; 
t[1] := FALSE; t{2] := FALSE; t[3] := FALSE;  t[4] := FALSE; 
t[5] := FALSE; t[6] := FALSE; 
fact[3] := computetrace(MAT , dim, 6, pos, t); 
pos(1] := Pl; pos[2] := U; t[1] := FALSE; [2] := FALSE; 
fact[4] := computetrace(MAT , dim, 2, pos, t); 
pos[1] := Pl; pos[2] := U; pos[3] := Pl; pos[4] := U; 
t[1] := FALSE; t[2] := FALSE; t[3] := FALSE; t[4] := FALSE; 
fact[5] := computetrace(MAT , dim, 4, pos, t); 
pos[1] := Q1; pos[2] := V; pos[3] := Ql; pos[4] := V; 
t[1] := FALSE; t[2] := FALSE; t[3] := FALSE; t[4] := FALSE; 
fact[6] := computetrace(MAT , dim, 4, pos, t); 
totaltrace := totaltrace + fact[{1] * fact[2] * fact(2] * fact[2); 
totaltrace := totaltrace + fact{3] * fact[4] * fact[4] * fact[4]; 
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totaltrace := totaltrace + 3 * fact[5] * fact(6] * fact[4] * fact[2]; 
totaltrace := totaltrace + 3 * fact[1] * fact[6] * fact[2]; 
totaltrace := totaltrace + 3 * fact[3] * fact{[5] * fact[4]; 
totaltrace := totaltrace + 4 * fact[1] * fact[3]; 


writeln(’The trace is: ’, totaltrace); 
END. 
By looking at the last 6 relevant lines and example 7.1.1, one can easily 
perceive the structure of this program. 

Someone who would like to compute this expectation for specific matrices 
P,, Q,, U and V, has to change the dimensions in line 2 and 3 and put the 
input matrices on a datafile. 

The module that generates this program has as an additional feature that it 
searches for equal factors. Instead of programming the computation of 20 
factors, only the 6 different ones are programmed. However, we do not pretend 
that the generated program is computationally efficient, since no time was 
spent on further refinements; we only use it to demonstrate a useful feature 


of formula manipulation by computer. 
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CHAPTER 8 


MATRIX ALGEBRA BY COMPUTER: MATRAL 


8.1 Introduction 


In chapter 2 an example of a very simple computer algebra system (VSCAS) has 
been presented. This system is capable of computing expressions assuming that 
the symbols in the expressions denote commutative variables. Matrices and 
vectors however do not have the commutative characteristic, so that VSCAS is 
not at all suited for performing theoretical computations containing matrix 
products. In this chapter a CAS — MATRAL - will be described that is able to 
perform matrix algebra. The matrix algebra it can perform is in its purest 
form: the symbols denote matrices and the matrices can not be filled with 
specific elements. Other computer algebra languages, like Maple or Reduce, 
have excellent options for performing matrix algebra on matrices of specified 
dimension and with specific non-numerical numbers. These CAS’s are capable of 
computing for example the inverse of such a matrix, but they lack the 
opportunity to handle matrix algebra easily, where no dimensions are given and 
without specific elements: matrix algebra where the symbols denote just 
matrices. In fields like econometrics matrix calculus is often used for 
theoretical derivations, where only some specific characteristic of a matrix, 
like symmetry, is known. Until now, these computations had to be done by hand. 

MATRAL is an effort to fill this gap in computer algebra. MATRAL knows 
among other things about a transposed matrix, an inverse matrix and the trace 
of a matrix expression. Also, matrices can be very easily defined to have a 
certain characteristic, such as idempotency or symmetry, and the user can 
specify that a matrix should be replaced by a matrix—expression, or vice 
versa. MATRAL will give as the output of an expression the simplified and 
expanded form of the input-expression; expanded using the given 
substitution-definitions and simplified using the characteristics (e.g. 
symmetry) of the matrices in the expression. MATRAL is also capable of writing 
an error-free Pascal-source for the computation of a freshly computed 
expression. 

In section 8.2 the capabilities of MATRAL will be described. In section 8.3 
some technical information will be given. Section 8.4 will tackle the problem 
of the computation of the variances of quadratic unbiased estimators of the 


error—component model with incomplete panels (Wansbeek and Kapteyn, 1989). 
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Section 8.5 concludes. | 


8.2 Description of the language of MATRAL 


First the user-interface of MATRAL will be discussed. Normally, each input 
line should start with a colon (:) followed by a. keyword. If this colon is 
missing, then MATRAL assumes that the input line contains a matrix—expression 
that should be computed. A summary of the keywords and capabilities of MATRAL 


will be given in this section. The following points of interest will be 


























discussed: 


a. Legal expressions. 

b. Characteristics for matrices. 

. Substitution of matrices. 

. Definition of matrices. 

. Simplifications of a matrix—expression. 


Traces and user-defined ordering. 


om Oo aa 


. Other keywords. 
Ad a. Legal expressions 


The context free grammar (see chapter 2 for more information), defining the 
set of all expressions which should be considered legal to our parser is given 


by: 


<expression> ::= <term> | <term> <operator> <expression> | 


<operator> <term> 





<term> = <factor> | <factor> ’| <term> <factor> | <number> 

<factor> = (<expression>) | <operand> 

<number>  ::= d | d <number> 

<operator> s:= + [-|/ 

<operand> = <matrix operand> | <non-matrix operand> 

<matrix operand> ::= v | v <number> 

<non-matrix operand> ::= w | w <number>, 
where v represents any character taken from the set [’A’ .. ’Z’]+[’] and w 
represents any character taken from the set [a .. °2’] - [#7]. The 
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apostrophe symbol ‘ denotes the transpose of a factor and a factor preceded by 
a slash indicates that this factor should be inverted. In MATRAL, a 
<non-matrix operand> is defined to have the commutative characteristic, 
whereas <matrix operand> is defined not to have the commutative 
characteristic. So variables starting with a capital character and the 
character i are assumed to be matrices (non—commutative variables), variables 
starting with a character from the rest of the alphabetical character set are 
assumed to be commutative variables. For example, the legal expression 
17a/bBC’ /DE/(2F +vG)i+HaQi represents 17 5 BC'D"E(2F +vG)i+aHQi. In 
MATRAL the matrix I is predefined to be an identity matrix of size n and i is 


predefined to be an n-vector of ones. 
Ad b. Characteristics for matrices 


In matrix algebra certain characteristics are commonly used. MATRAL knows the 
function of the characteristics noted below. For example, the user can define 
a matrix M to be idempotent. Then each time that an expression containing MM 
is found, MATRAL will simplify this into M. The table below will show the 
known characteristics in the current version of MATRAL and an example how to 
learn MATRAL that a specific matrix has a certain characteristic. 


Characteristic Input form for Matral 
1. Identity identity <M> [for <M1>] 
2. Symmetric :[+|-] symmetric <M> 
3. Idempotent :[+|-] idempotent <M> 
4. Vector of ones :onevector <M> dimension <n> 
5. Projector on a certain matrix. :[+|-] projector <M1> on <M2> 
6. Zero :zero <M> 
7. Diagonal ‘diagonal <M> 


Note: | means or, [] means optional, <> means obligatory, <M> means a 
matrix operand, <n> means a non—matrix operand. A + before a keyword means: 
add characteristic to previously defined characteristics. A —- before a keyword 


means: remove from previously defined characteristics. 
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Example 8.2.1 


Input for MATRAL Output from MATRAL 

sidentity I for Q I identity matrix for og? 

:symmetric Q Q is symmetric 

:symmetric P P is symmetric 

:tidempotent P P is symmetric, idempotent 

:+projector P on Q P is symmetric, idempotent, 
: projects Q to 0 

(I+P+Q)(I+P+Q)’ I+3P+20+00 


So given the characteristics above, MATRAL expands and _ simplifies 
(I1+P+Q)(I+P+Q)' into 14+3P+20+Q0. 


Ad c. Substitution of matrices 
In MATRAL a matrix can be defined to be a matrix expression. This means that 
each time MATRAL encounters this specific matrix in an expression, it will be 
replaced by the defined matrix expression. The syntax in MATRAL is: 
:substitute <M> = <expression>, with <expression> any legal expression. 
Example 8.2.2 
:substitute Q = I-X(X’X)"X’ 
MATRAL will now replace Q in an expression with [-X(X nS ae, each time Q is 
encountered. So [+3P+20+Q0 will yield 41+3P-3X(X'X yx '. Matral does not know 
that for this specific Q it holds true that QQ=Q, but it just simplifies the 
terms and collects identical terms. 
Ad d. Definition of matrices 
In MATRAL a certain group of matrices can be replaced by another matrix, which 
is the complement of point c. Each time MATRAL encounters this group of 
matrices it will be replaced by a single matrix. The syntax in MATRAL is: 


:[+|-]define <M> = <single—expression>, 


with <single-expression> such that the expanded form yields an expression 
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without additive terms. E.g., 


:define X2 = X(X+Y)"'Z is allowed whereas 
:define X2 = X+Y is not. 


Example 8.2.3 

:define X1 = X(X'X) 1X’ 
MATRAL will now replace X(X ass Teg by X1 each time it is encountered. So 
41+3P-3X(X'X)"X’ will yield 41+3P-3X1. 


Ad e. Simplifications of matrix — expressions 


Obviously, expressions will be simplified using the characteristics and 
definitions given to certain matrices. Also equal terms will be found and 
added, zero terms will be deleted. For a discussion on how to find equal 
terms, see chapter 2. More simplifications can be done: MATRAL detects and 
simplifies certain expressions involving inverses of matrix—expressions. For 
example: A"'A , AB(AB)™*, A(A+B)"'+B(A+B)"! will all be simplified to the 
unity matrix I. Also a more complicated term like, 


BDB'A"B(B'A"B + D*) "BA? + B(B'A™B + D?)"'B A 


will be simplified to BDB’A™. In the next section the algorithms used by 
MATRAL to simplify such expressions will be given. 


Ad f. Traces and user-defined ordering 


MATRAL is capable to compute the trace of a matrix—expression. The syntax 


activate or deactivate computing according to the trace rules, is 
:[+|-]trace 
The following rules are used to simplify trace expressions. 


1. tr(AB)=tr(BA) 
2. tr(A)=tr(A’) 


A unique ordering is used (see also chapter 7) in order to detect equal terms. 
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By default, output will be generated using the alphabetical or ASCII-ordering. 
The resulting output may not please the user: for example if the resulting 
expression has to .be transformed into a computer program that handles real 
data, then the most efficient expression should be given as output, in order 
to minimize the computing effort. 

For example, let Y be a txh matrix and X an Axh matrix, with t << h. 
Actually computing the matrix products in tr(XX’Y’Y) from left to right would 
cause O(n’) multiplications. However computing the matrix—products in the 
equal expression tr(YXX'Y’) from left to right would cause only O(h?) 
multiplications. So, even though tr(XX’Y'Y) = tr(YXX’Y’), the output-form 
matters when actually computing the resulting expressions with real data. 

Because it is desirable that MATRAL produces a more efficient output form 
than the ASCII-ordering, this ordering can be overruled by teaching MATRAL a 
new alphabet. This alphabet consists of <matrix operands> which should precede 


the ASClII-ordering, thereby also indicating the transpose-sign. The syntax is: 
salphabet <M1>[<M2>[...]] 


Example 8.2.4 
salphabet 7’Y 
:+trace 
XX 4+ XX'V'Y 
The answer will be tr(i’XX’i)+tr(YXX'Y’). 


Ad g. Other keywords 


A few keywords have not been discussed. The list below gives them with an 


explanation of their function. 


:clear <variable> clear a definition of a variable 
sallclear clear all variables 

:compute <expression> compute the expression 

:‘recompute recompute the last result 

pascal write a Pascal program 

sstore <M> store the last computed result in <M> 
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8.3 Technical description of MATRAL and some algorithms 


In MATRAL the expressions are stored in expanded form. Since MATRAL is also 
capable of holding expressions of both the commutative type and the 
non-commutative type, this is done in a data structure as_ indicated 


schematically in diagram 8.3.1. 


Diagram 8.3.1 

sca|comm. |non— sca|comm. ]non- 
Lp [ise| "Yen, [-}-* ~- — [ise] [corn 
In this diagram scalar means the scalar part, comm. means the commutative 
part, and non—comm. means the non—commutative part of the expression. The 
arrows (pointers) mean an addition and each block represents the 
multiplication of scalar, the commutative part and the non-commutative part. 
Since the main task of MATRAL is storing and manipulating matrix—expressions, 
and the scalar part and the commutative part have already been discussed in 
chapter 2, only the non-commutative part will be discussed below. The 
commutative part has to be able to store matrix—expressions, possible 
consisting of the transpose sign and/or inverse-sign, for example the 
expression 2AC’(Y'X+aR)'Z +Q. 

With only our concentration on the non-commutative part, this expression 
can be stored schematically as shown in diagram 8.3.2. The arrows pointing 
horizontally indicate additive terms, the arrows pointing vertically indicate 
a multiplication or an inversion, depending on the contents of the first 
element in the record it points to. 

From now on we will use the word term for each record with possibly a 
pointer to an additive term, where the pointers to an additive term will be 
ignored; expression for the combination of one or more terms, so the pointers 
to additive terms are not ignored; and factor for a record that may contain an 


element or a _ pointer to an _  inverse-expression. In diagram 8.3.2 


2) | |= is an example of a term, elaml [clam | is an 
example of an expression, and Baral and Ml | are both examples 


of a factor. 
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Diagram 8.3.2 


a 








Car ie a a) 


BREAD 










BREAD 


In Pascal the definition of a data structure as indicated in diagram 8.3.1 and 
8.3.2 could be like this: 


type 
matfactorptr = “matfactor; 
matfactor = record 
next: matfactorptr; 
Case divide_symbol: boolean of 
TRUE : (inverse: mat_termptr); {corresponds with /} 
FALSE: (main: mat_ element); {corresponds with +} 
end; 
mattermptr = “matterm; 
matterm = record 
scalar : scalartp; {not discussed} 
comm : commutative_factorptr; {not discussed} 
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noncomm: matfactorptr; 
next : mattermptr; 
end; 
matelementptr = “matelement; 
matelement = record 
mat: word; {matrix-symbols are converted to numbers} 
transpose: boolean; 
end; 


Below some algorithms will be given for simplifying expressions like 
CABB™*(I+A)"+C(I+A)* to the expression C. This type of simplification will 
be called inverse simplification. The algorithm for finding all expressions 


that can be simplified by inverse simplification, is divided in three parts: 


1. search-for—inverse—simplification—in—all—terms 
2. search-for—single-term—inverse—simplification 
3. search-for—multiple-terms—inverse-simplification 


We will first give the algorithms and then show how these algorithm work by 


giving an example. 
Algorithm 8.3.1 search—for —inverse—simplification—in-all-terms 


1. [Initialize] e is a pointer to the expression to be evaluated. Let head 
point to e. 

2. [Loop] If e = nil then go to step 9 otherwise let t <— e. 

3. [Get inverse] Let f be the first factor of t which points to an 
inverse—expression, i.e. in diagram 8.2: where divide-symbol = ’/’. 

4. [Inverse—-expression exists?] If f = nil, ie. in this term ¢ no inverse— 
expression is available, then go to step 8. 

5. [Test] Let v be the pointer to the inverse-expression belonging to f, 
ie. f*.inverse. First it will be checked if this expression can be 
simplified itself by inverse simplification. This will be done by a 
recursive call to algorithm 8.3.1, where in this case e — the expression 
to be evaluated - is equal to v. After completion, the algorithm will 
resume at step 6. 

6. [Single or multiple] If v is a single term, i.e. consists of exactly one 


term, then call algorithm 8.3.2 otherwise call algorithm 8.3.3. 
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[Next inverse] Let f be the next factor of ¢ which points to an 
inverse—expression, i.e. in diagram 8.2: where divide-symbol = ’/’, and 
go to step 4. 


[Next term] e <— e*.next and go to step 2. 


9. [Finish] Restore initial value of e: e < head. 


Let copy(t, i, j) be a copy of j factors starting at factor number # of term 


t. Let factor_len(t) be the number of factors in term ¢, counting a pointer to 


an inverse expression as one, and let term_len(e) be the number of terms in 


expression e, counting only the main additive terms. 


Algorithm 8.3.2 search—for—single-term—inverse—simplification 


[Initialize] v is the pointer to the inverse-expression, which is a 
single term. Compute c3 <— factor_len(t), c, <— cz — factor_len(f), 


Cz <— factor_len(v). 


2. [Identity?] If v represents an identity matrix, then go to step 9. 


3. [Enough factors in front?] If cz > c, then go to step 6. 


Now the following steps will search for a matrix in front of v, for 
example: AB( AB)’. 
p <— copy(t, cy-cyt1, ¢2). p now contains exactly as many factors as v. 


5. [Equal?] If p is equal to v, except for a constant and/or a commutative 


term, then go to step 9 otherwise go to step 6. 
[Enough factors in back?] Now the following steps will look for a matrix 
behind expression v, for example: (AB) AB. 


If cp > (¢3-€,-1) then go to step 10. 


7. p < copy(t, c,+1, c2). p now contains exactly as many factors as v. 


10. 


{Equal?] If p is equal to v, except for a constant and/or a commutative 
term, then go to step 9 otherwise go to step 10. 

{[Simplify] Term ¢ can be now be simplified, so simplify term ¢ (not 
explained how). 

{Exit] Finish scanning this expression v, and return to the calling 
algorithm. 
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Algorithm 8.3.3 search—for—multiple-terms—inverse—simplification 

The first part of this algorithm will look for pieces of the inverse vin 
front of v. If a ”candidate” has been found, the expression e is searched for 
the other pieces of the inverse. For example, if the expression e equals 
BC(A+C)'A+BC(A+C)'C then v equals (A+cy". The candidate that might be used 
to simplify the expression, is C(A+C)". The rest of the expression will be 
scanned for the missing pieces (in this case: BA(A+C)"A). If not all missing 
pieces are found, the algorithm will look for pieces of the inverse v behind 
v. In this example, the candidate would be (A+C) "A. Now the rest of the 
expression will be scanned for the missing pieces (in this case: BC(A+C)'C). 
Since the second part of this algorithm is similar to the first part, only the 


first part will be given. The pieces will be stored in array pce. 


1. [Initialize] n, <— term_len(v), / < 1, 
¢, <~ number of factors in t before factor f was encountered. 

2. [Possible?] If c, = 0 go to step 14. 

3. [Loop] If | > n, then go to step 14. 

4. [Initialize] Let v, be the [-th term of expression v. If v, represents 
an identity matrix then d<O and x <— v,, otherwise d <— factor_len(v) 
and x <— copy(t, ¢,-d + 1, d). 

5. [Equal?] If x # v, then go to step 13. 

6. [Yes, equal] 9g, <- copy(t, 1, c¢-d), go < v, peli] < e,. If /=1 then 
k«~—2 otherwise k<—1. We have found what could be a piece of the inverse. 
Now we will scan expression e for the remaining pieces. 

7. [Initialize] e,<e, ok < false. 

8. [Compute] Compute e, < g, * v, * 92. 

9. [Look for e, in e,} If e,; is already stored in pc[m], m=1,..,k-1 or in 
pe[t], then go to step 10. If e, equals e, then go to step 11. 

10. {Loop} e, < e,".next. If e, # nil go to step 9 otherwise go to step 13. 

11. [Part found!] ok <— true, pce[k] < e,, k << k + 1. If k=l then ke—k+1. 

If k < n, then go to step 7. 

12. [Total found!] Now we have found all the necessary n, pieces. We can 
simplify expression ¢ by eliminating from the expression in pe{l] v*g2 
and eliminating all other pieces. Go to step 13. 

13. [Loop back] | <~ | + 1. Go to step 3. 

14. [Exit if done] If ok = true then exit: simplification has been done, 


otherwise look for pieces behind v (not discussed). 
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Example 8.3.1 
This example will simplify the expression CABB (I +A) SCI tay using the 
algorithms 8.3.1-8.3.3. Diagram 8.3.3 shows how this expression is stored in 


memory. 


Diagram 8.3.3 Representation of CABB™ TUI+A yc +A)” 


eo cart 7 


= 
Hla. ,, Ell, 


DBE. ial it i 
"76 
COE 








Lei, il | 4 |e 
Bio B 


In this diagram the blocks of information are denoted by B;. These blocks will 
be used to indicate which information is erased from memory or to indicate the 
start of a pointer variable. A chain of B,’s is used to indicate the contents 
of a variable, for example, By—>B,o—B,;>B,, denotes I+A. We will now 
simplify the expression above by giving the results’ of the relevant steps. 


Algorithm 8.3.1 
Step Result 
1,2,3,4 e points to B,, t points to B,, f points to Bs, v points to Bg. 
6 v is a single term —> 
Algorithm 8.3.2 
1 v points to Bg, ¢<—3, cz<—1. 
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4,5 pis a copy of B, and equals v. 


9 


t will be simplified: the blocks By, Bs, Bs, Bz will be 
erased. Return to algorithm 8.3.1. 


f points to By. 


v points to Bg. 


v is consists of more than one term —> 
Algorithm 8.3.3 


1 
4 


aonrn vu 


10 


10 
13 


10 


11 
12 


Nye—2, lel, Cy<—2. 

v, points to Byg; since it represents an identity matrix, 

d <— 0, and x is a copy of Byp. 

xX <— Vy. 

9, is a copy of B,—>B3, g, points to v, pc[1]<e, k<—l. 

€, points to B,, ok <— false, k <— 2. 

C2 <— gy * V2 * Go = By—>B3—>Byz By By By Bp, 

which represents CAA(I +A)”. 

e, is already stored in pe[1]. 

€,; now points to B,3. 

€2 is not equal to e,. 

e, < nil. 

l<—2, go to step 3. 

v2 points to By, d <1, x is a copy of B3. 

x <— %. 

9, is a copy of B,, gz points to v, pc[2] points to B,, 
ke—1. 

e, points to B,, ok < false. 

€2 <— 91 * Vy * Go = By—>By—>By—Byy By B2, 

which represents C(I+A)”. 

€, is already stored in pc[2], so go to step 10. 

e, points to By. 

€, equals e,. Go to step 11. 

ok«—true, pc{1] points to B,, k <— 3, k > n, (3 > 2) 

All necessary pieces have been found: the expression e can 
be simplified by eliminating the blocks B3, Bg, By, By, 
By, By, from pe[2] and by eliminating pc[1]..So expression 
e is simplified to B,—>B,. Return to algorithm 8.3.1. 


f< nil, e — nil. 


e is simplified to B,—>B, and represents the matrix—expression C. 
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8.4 Examples 


As a test-case for MATRAL the variance of quadratic unbiased estimators of the 
error—component model (ECM) with incomplete panels will be computed. Wansbeek 


and Kapteyn (1989) describe the following regression model: 
Yt = XntVet XB + Une, (8.4.1) 


where h (hk = 1,..,4) denotes households and ¢ (t=1,..,7) years; x,, is a 
k-vector of explanatory variables; @ is a k-vector of parameters; the 
error-term w,,; are iid. random variables with mean zero, and independent of 
the Xpz- 

In the fixed—effects (FE) model o, and y, are fixed constants. In the 
random—effects (RE) model o;, and y; are i.id. random variables with mean 
zero, mutually independent and independent of the x,, and tng. 

Let N, (N; < H) be the number of observed households in year ¢ and N=E,N,. 
Let D, be the N,xH identity matrix from which rows corresponding to households 


not observed in year t have been omitted and consider 


D, Dyiy 
Z = (Z, %) =| : % (8.4.2) 


The matrix Z gives the dummy-—variable structure for the incomplete-data model: 


model (8.4.1) can be rewritten in vector form: 

y= Za+ Zoey + XB +4. (8.4.3) 
Let 

Ay = 232Z,, Ap = 2922, A=Z3Z), (8.4.4) 
then Ay is a diagonal HxH matrix with the h-th element indicating the number 
of years for which the h-th household has been observed, Ay is a diagonal 
(TxT) matrix with the t-th element indicating the number of observations in 


year t and A is a TxH matrix of zeros and ones indicating the absence or 


presence of a household in a certain year. Further define, 
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2 = 2-24, M (8.4.5) 
Q = Ay-AA; A’, (8.4.6) 
Pm ly = ZAG = 207’, (8.4.7) 


where Q is a generalized inverse of Q. 

Wansbeek and Kapteyn (1989) show that P is the projection matrix onto the 
null-space of Z. They also show how to transform the FE model with incomplete 
panels in order to estimate the regression coefficients: the OLS—estimation of 
the regression coefficients in model (8.4.1) is equal to OLS—estimation of the 


regression coefficients in model: 
Py = PXB + t. (8.4.8) 


Let o? = var(upz), 07 = var(o,), 0% = var(y;), then the covariance of the 


composite error term Ep, = Up, + O; + ¥¢ is 

Q = oly + 072,Z;, + 032.2). (8.4.9) 
A way to efficiently estimate the ECM is to estimate the variance components 
o*, o?, and o% by a quadratic unbiased estimation method (QUE) and then 


estimate the regression coefficients by GLS with these estimates inserted in 


92. Now we can have two situations: 


1. The matrix X of explanatory variables does not contain a vector of ones. 


2. The matrix X of explanatory variables contains a vector of ones. 
8.4.1 The matrix X of explanatory variables does not contain a vector of ones 


Let e be the N—vector of FE residuals, i.e. e=y—-Xb, with b the FE OLS—estimate 
of B. So e=y-X(X'PX)"X'’Py = My = Me. Now define 


Oy = CZ,Ar Ze, (8.4.10) 
Gr = e'Z,Aj Zie, (8.4.11) 
ky = tr((X’PX) 'X'Z,Az'Z}X), (8.4.12) 
kp = tr((X'PX)"X'Z,Aj, ZX). (8.4.13) 


Then it can be shown — for example by MATRAL - that 
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E(qu) = (T+ky) 0? + To} + No, (8.4.14) 
E(qr) = (H+kp) 0? + NoZ + Ho}, (8.4.15) 
E(e'Pe) = (N-T-H+1) 0. (8.4.16) 


These three equations in the three unknowns o?, of and o% can be solved, which 


lead to the following estimations: 


6? = e’Pe / (N-T-H+1), (8.4.17) 
82 = (e'Pe(NH-TH +Nkp—Nky)/(N-T-H+1)+Hqy—Nar) (HT -N’), (8.4.18) 
82 = (e’ Pe(NT -TH +Nky—Tkr)/(N-T—H+1)+T¢r—Nay) /(HT-N’). (8.4.19) 


Let @ denote any of the variance components o?, 0? and o%. It is clear that 


quadratic unbiased estimators can be written as 

6 = e'We = e'M'WMe, (8.4.20) 
(Wansbeek and Kapteyn (1989, p.356) erroneously state 6 = e'MWM'e) 
with 

W = aP4bZ,Aq'Zi+cZ,Az Z}. (8.4.21) 


The constants a, 6 and ¢ can be chosen such that 6 is unbiased, i.e. the 
relevant formulas from (8.4.17)-(8.4.19) hold. If we assume that €,;, ©, and 
¥, are normally distributed, we are able to compute the variance of 6: 


according to the multivariate normal theory, 
var(6) = 2 E(tr(e’M’WMee'M'WMe)) = 2 tr((M’WM2)’). (8.4.22) 


Elaborating this formula by hand is highly tedious. MATRAL is very well suited 
for performing manipulations to find the answer to formula (8.4.22). The input 
for MATRAL as given in table 8.4.1 will let MATRAL elaborate this formula. For 
readability reasons some symbols in the table are used which are in fact not 
allowed in MATRAL. 
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Table 8.4.1 Input for MATRAL 


1.> :projector P on Q 

2.> :+projector P on 2, 

3.> :+projector P on Z, 

4.> :+symmetric P 

5.> :+idempotent P 

6.> :+symmetric Q 

7.> :+symmetric S 

8.> identity I 

9.> :+define I = SX'PX 

10. > :substitute W = aP+Q 
11. > :substitute 2 = 0?I+03?Z,Z\1+03Z.Z, 
12.> :substitute G = I 

13.> :substitute M = I-XSX'P 
14.> :compute M'GWGM2Q 


16.>:recompute 
17.>:store in T, 
18.>:+trace 
19.>:compute T,7T, 
20.>:+identity I, for 2, 
21.>:+identity I, for Z, 
22.>:symmetric Ay 
23.>:symmetric Ar 
24.>:+define Ay=Z;Z, 
25.>:+define Ap=Z322 
26.>:recompute 
27.>:define A=Z32, 
28.>:alphabet X'Z,AA 
29.>:recompute 

















15.> :define Q = b2}/(2;2,)Z,+¢Z/(Z2Z2)Z2 
Using this input MATRAL computed formula (8.4.23). 


var(6) = 20703c? tr(X’Z, Z; X S) (8.4.23) 
+4o%o3be tr(X’Z, A Ay'Z} X S) 
4+2070%¢? tr(X’ Z, Ap'Z; X S) 
+ o%0%c? tr(X’ Z, Ap'Z, X SX’ Z, Az'Z, X S) 
+2070%be tr(X' Z, Ap'Z, X S X' Z, Ay'Z; X S) 
+207o}c? tr(X’ Z, Ap'A A’ Ap'Z; X S) 
+40%0?be tr(X’ Z, Ap'A Zi X S) 
+407obe tr(X Z, Ap A Ay'Z;'X S) 
+20702b? tr(X’ Z, Z, X S) 
+20703b? tr(X’ Z, Ay’ A’ A Ay'Zj X S) 
+20%07b? tr(X’ Z, Ajy' Zi X S) 
+ 0%07b? tr(X’ Z, Ay Z; X SX’ Z, Ay'Z; X S) 
+ (20}03b? + 4o?o03be+2020%c? ) tr(A A’) 
+ (4070?bc+2070%?c?) tr(A A’ Az’) 
+ ofo3c? tr(A A’ Az'A A’ Az’) 
+2o30%be tr(A Ay A’ Az) 
+(20203b? +4c03bc) tr(A Ay’ A’) 
+0303b? tr(A Ay A’ A Aj’ A’) 
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+20}03be tr(A Aq’ A’ Ar) 
4+20%02be tr(A Aq A’ Az’) 
4+207a%c? tr( Ar) 

+o30%c? tr(Ar Ar) 
+20703b? tr( Ay) 

+03202b? tr(Ay Ay) 
—a?o7o? tr(I) 

+o0707b? tr(I,) 

+0702c? tr(I2) 

+a?0?o? tr(P), 


where S=(X’PX)"1. This term can be simplified by hand, by noting that tr(I)=N, 
tr(1,)=H, tr(I,)=T, tr(P)=N-H-T+1, tr(Ay) = tr(Ar) = N. 


8.4.2 The matrix X of explanatory variables contains a vector of ones 


The derivation of the model where the matrix X of explanatory variables 
contains a vector of ones, is similar as under 1. Formula (8.4.10) and 


(8.4.11) change into 


PEAp oe, (8.4.10’) 
POA es, (8.4.11') 


qH 
qr 


with f = Exe and Ey = I - aii’. Now it can be shown - for example by MATRAL - 
that 


E(Gu) = E(u) — (1+éyX(X'PX)"')o? + inZ,Zjino? + inZ2Zaino} — (8.4-14’) 
E(Gr) = E(qr) — (1+iyX(X'PX)*)o? + inZ,Zjiyo? + inZaZaivo} — (8.4.15") 


Let 6 denote any of the variance components 0, o? and o3. Analogously to the 


first model, the quadratic unbiased estimators of this model can be written as 
6 = fWf = ©'EyM'WEyMe, (8.4.20') 
with 


W = aP+bZ,Ay'Zi+cZ,Az 23. (8.4.21) 
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Of course, the constants a, b and ¢ are not the same as in (8.4.21). They can 
be chosen such that (8.4.20’) is unbiased. Once again we are able to compute 
the variance of 6: 


var(9) = 2E (tr(e’EyM 'WEyMee'EyM 'WEyMe)) = 2 tr((M' EyWEyMQ)*)(8.4.22') 


By adding and modifying the following lines to the previously stated input for 
MATRAL, 


12. >:substitute G=I-it' /n 
12a.>:onevector i dimension n 
12b.>:+define i=Z,/(ZjZ,)Zii 
12c.>:+define i=Z,/(Z23Z2)Z3t 
28. >:alphabet #’X'Z3,AAr 


formula (8.4.22’) can be computed. The result is: 


var(9) = (8.4.23) + (8.4.23") 
( -0707b? - 20707bc — 0702?) 
+ ( —20707b? — 40207 be —2070%c?)/n i’ XS Xi 
+(0707b +20707bc+070%c?)/n2 vw XS X iv XSK i 
+ (20703? + 407 03bc + 20703c?)/n? i’ XS X' it’ Z, Z' i 
+ (20202b? + 40702bc + 2020%c?)/n2 i’ X SX ii Z, Zi 
+(-4070%be-407a3c?)/n i’ X SX’ Z, Z,' i 
+(-20702be-2070%¢2)/n i X SX’ Z, Ap Zz, X SX i 
+(-40%0?be— 40703c?)/n i’ XS X' Z, Ap AZ i 
+ (-40702b? —40702bc)/n it XS X' Z, 2, i 
+ (—40703b? ~ 40703be)/n i’ XS X' Z, Ay’ A’ Zy' i 
+(—20702b? -2070%be)/n i’ X S X' Z, Ay’ Z,' XS Xi 
+ (-20703b? - 407 03be — 20703c?)/n i Zy Z,' i 
+ (0}03b? + 20303bc+ o303c?)/n? i Z, Zo i i’ Z, Zo’ i 
+ (2070%b? + 4070%bc + 2a70%c?)/n? i Z, Z,' i i Z, Z,' i 
+( -40f0%b? - 80203be ~40703c?)/n i’ Z, A Z,' i 
+ ( —20303b? -20%0%bc)/n i’ Z, A Ay’ A’ Z,' i 
+ ( —20%03bc - 20303c?)/n t’ Z, Ap Zy' i 
+(-0703be-0703c?)/n i’ Z, Ap AZ,’ i 
+ ( -2020?be — 2020c?)/n i’ Z, A’ Ap’ AZ,’ i 
+ ( ~20303b? - 20302bc)/n # Z, Ay Zy' i 
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+ ( -20703b? — 307 a?bc—07a%c?)Jn i Z, Zt 


+ (0207b? + 2020?bc + o207c?)/n? i Z, Zy) it Z, Zy' i. 


It would be rather bothersome to actually program the results (8.4.23) and 
(8.4.23’). Wansbeek and Kapteyn state that they have actually done this. If 
MATRAL would have existed at the time they programmed these formulas, much | 
time could have been saved: MATRAL can write Pascal source code for matrix 


formulas as given above. 


8.5 Conclusion 











In this chapter we have seen a computer matrix algebra system 5 MATRAL — which 
can simplify matrix expressions based on a priori known characteristics of 
matrices. MATRAL knows” about certain often used characteristics, like 
symmetry or idempotency, and simplifies expressions accordingly. This system 
should be seen as a tool for theorists to compute formulas, given some basic 
information about the various matrices in the expressions to compute. The 
computation done by MATRAL is just a combination of expanding and simplifying 
a given matrix expression. It could also be used for practical purposes: for 
the programming of a rather annoying matrix formula. By just entering the 
formula and pushing one more button, a flawless computer program will be 
given, which can compute the formula for real data. 

MATRAL is still rather a primitive system: many options are missing. For 
example, the vec-operator and the Kronecker product are still unknown to 
MATRAL. These and other options will be added in future versions of MATRAL. 
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CHAPTER 9 


CONCLUSION 


The aim of this study was to integrate computer algebra in econometrics and 
statistics. It is shown in chapters 3, 4 and 6 that computer algebra is a 
useful instrument for characterizing the identification situation in 
simultaneous equation models with measurement error, restricted Factor 
Analysis models and general LISREL models. Without the use of computer 
algebra, this characterization is almost impossible to handle. In chapter 5 
another computer algebra program — MINPARAM -is used to find automatically a 
parameterization of the reduced form of simultaneous equation models. In 
chapter 7 a theoretical statistical problem of the expectation of higher- 
degree functions of normally distributed matrices is solved by a computer 
algebra program. A freshly computed formula can be ’programmed’ automatically 
in Pascal source code so the expectation involved can be computed easily for 
real data. In chapter 8 a computer algebra system —- MATRAL -— is given that 
handles matrix—expressions. We show that this system can be a useful tool for 
the theoretical econometrician in solving certain analytical problems. 

Of course computer algebra applications are not limited to these subjects: 
anything that can be done by algebra with paper and pencil can in principle 
also be done by computer algebra. Computer algebra per se is not an 
econometric or statistical subject: it belongs to computer science. For 
econometricians computer algebra should be seen as a tool, that can be used to 
solve analytical problems which are complicated or tiresome to perform by 
hand, but rather straightforward in the repeated application of a limited set 


of rules. 
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Computeralgebra programma’s zijn in staat om symbolische wiskunde te 
bedrijven; ze kunnen niet alleen numerieke expressies uitrekenen zonder 
afrondfouten, maar ook formules manipuleren die nog geen numerieke waarden 
hebben; elk algebraisch object kan exact in het geheugen van de computer 
opgeslagen worden. 

Computeralgebra wordt tot op dit moment vrijwel niet gebruikt om 
econometrische problemen op te lossen, hoewel het een bruikbaar instrument kan 
zijn voor zowel de toegepaste als de theoretische econometrist. Een toegepaste 
econometrist kan bijvoorbeeld een computeralgebra systeem (CAS) de afgeleide 
laten berekenen van een bepaalde doelfunctie, en het resultaat automatisch en 
foutloos laten programmeren. Een theoretische econometrist kan een CAS 
gebruiken om ingewikkelde of omvangrijke manipulaties uit te voeren op 
analytische expressies. 

Dit proefschrift behandelt de toepassing van computeralgebra in de 
econometrie en statistiek. Met andere woorden, er wordt bekeken of 
computeralgebra een handig instrument is om een verscheidenheid van 
econometrische of statistische problemen op te lossen. Er wordt nauwelijks 
aandacht geschonken aan computer algebra systemen: alle programma’s die 
behandeld worden, zijn geschreven met behulp van de programmeertaal Pascal. 
Een voordeel van het gebruik van Pascal is dat beter inzicht verkregen wordt 
in de beperkingen van computeralgebra en dat er geen beperkingen zijn in de 
manier van implementatie: elke manipulatie en elk algorithme kan worden 
geprogrammeerd; in computeralgebra systemen kan een bepaald algorithme niet 
aanwezig zijn, of een bepaalde bewerking ongedefinieerd zijn. Een voordeel van 
het gebruik van een CAS is dat een gebruiker niet hoeft te weten hoe symbolen 
worden opgeslagen in het geheugen van de computer en dat er 
voorgeprogrammeerde algorithmes aanwezig zijn. 

Hoofdstuk 2 bespreekt hoe computeralgebra programma’s geprogrammeerd kunnen 
worden in een hogere programmeertaal, zoals Pascal. Hiertoe worden alle 
ingredienten gegeven om een simpel computeralgebra systeem te bouwen; de 
appendix bevat ook daadwerkelijk de broncode van zo’n CAS. 

Hoofdstuk 3, 4 en 6 behandelen het identificatie-probleem van verschillende 
modellen. In al deze hoofdstukken wordt er vanuit gegaan dat de waarnemingen 
normaal, identiek en onafhankelijk verdeeld zijn. Het identificatie-probleem 
komt dan overeen met de vraag of de parameterwaarden uniek bepaald kunnen 


worden uit de eerste— en tweede-orde momentvergelijkingen. Over het algemeen 
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worden rangcondities gegeven voor de identificatie van de parameters van een 
structureel model, waarbij niet wordt gezegd hoe de toegepaste econometrist de 
betreffende rangconditie moet evalueren. In hoofdstuk 3 wordt een 
computerprogramma (ERA) gepresenteerd dat de rang en nulruimte kan berekenen 
voor een verscheidenheid van modellen. 

Hsiao (1983) en Geraci (1976) hebben rangcondities gepresenteerd die het 
identificatie-probleem van simultane _vergelijkingen met meetfouten 
karakteriseren. In hoofdstuk 3 wordt aangetoond dat evaluatie van deze 
rangcondities tot de verkeerde conclusie omtrent de identificatie kan leiden. 
Een alternatieve rangconditie wordt afgeleid en een computerprogramma (ERASIM) 
wordt gepresenteerd dat de restricties waar de parameters aan onderhevig zijn, 
gebruikt om de identificatie situatie van de individuele parameters te 
karakteriseren. 

In hoofdstuk 4 wordt de identificatie van gerestricteerde factor analyse 
modellen beschouwd. Er wordt ook hier een algemene rangconditie afgeleid en 
een computeralgebra programma gepresenteerd (ERAFAC) dat automatisch de 
Jacobiaan kan construeren en evalueren om de identificatie van de individuele 
parameters te karakteriseren. 

Bekker en Dijkstra (1990) bespreken het exacte aantal restricties dat op de 
gereduceerde vorm van een simultaan model met exclusie restricties rust. In 
hoofdstuk 5 worden een algorithme, grotendeels door Bekker afgeleid, en een 
computer programma (MINPARAM)  gepresenteerd die deze restricties 
karakteriseren en een minimale parameterisatie van de gereduceerde vorm 
leveren, d.w.z. een parameterisatie met zo min mogelijk parameters waarbij de 
restricties op de gereduceerde vorm in beschouwing worden genomen. 

Een minimale parameterisatie van een matrix zoals Br wordt ook in 
hoofdstuk 6 gebruikt. Dit hoofdstuk behandelt het identificatie-probleem in 
LISREL modellen (Jéreskog, 1984). Een aantal rangcondities worden afgeleid en 
ook een aantal computeralgebra programma’s worden gepresenteerd die deze 
rangcondities kunnen construeren en evalueren. Het eerste deel van dit 
hoofdstuk behandelt het identificatieprobleem van specifieke submodellen van 
het algemene LISREL model, terwijl het laatste deel van dit hoofdstuk het 
identificatie-probleem van LISREL modellen behandelt met algemene lineaire 
restricties op de parametermatrices. 

In hoofdstuk 7 wordt het theoretische probleem, van de verwachting van 
functies van normaal verdeelde matrices besproken. Over dit onderwerp bestaat 
een grote hoeveelheid recente literatuur, zie bijv. Neudecker en Wansbeek 


(1987). Er wordt aangetoond dat de verwachting van deze functies algebraisch 
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geévalueerd kunnen worden door een methode die repeated conditioning wordt 
genoemd. Er wordt een computer programma gepresenteerd dat een grote 
verscheidenheid van dit soort problemen kan oplossen en tevens in staat is om 
de zojuist berekende formule in Pascal source code om te zetten. 

Hoofstuk 8 introduceert een soort kladblok dat in staat is om matrix 
algebra uit te voeren, in de zin dat symbolen matrices voorstellen en de 
matrices niet gevuld kunnen worden met specifieke elementen. Matrices kunnen 
bepaalde karakteristieken hebben, zoals symmetrie, of een matrix kan vervangen 
worden door een matrix expressie. Dit kladblok — MATRAL - kan gebruikt worden 
om een diversiteit van symbolische matrix—expressies te berekenen. Hier wordt 
het toegepast om de variantie van kwadratische zuivere schatters van het 
error—component model met incomplete panels (Wansbeek en Kapteyn, 1989) te 
berekenen. 

Dit onderzoek toont aan dat het gebruik van computeralgebra een handig 
hulpmiddel is voor het oplossen van enkele econometrische problemen. Het nut 
van computeralgebra blijft niet beperkt tot de problemen die hier behandeld 
worden: in principe kunnen alle wiskundige manipulaties die met pen en papier 
uitgevoerd worden ook m.b.v. computeralgebra gedaan worden. Computeralgebra op 
zich zelf is niet een econometrisch onderwerp; het behoort tot de informatica. 
Voor econometristen moet het als een hulpmiddel gezien worden dat gebruikt kan 
worden om analytische problemen op te lossen die ingewikkeld zijn, of 
vervelend zijn om met pen en papier aan te pakken. 
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